{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import sklearn\n",
    "import numpy as np\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn import preprocessing\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tteofili/.local/lib/python3.6/site-packages/IPython/core/interactiveshell.py:3063: DtypeWarning: Columns (9) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  interactivity=interactivity, compiler=compiler, result=result)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LoanNr_ChkDgt</th>\n",
       "      <th>Name</th>\n",
       "      <th>City</th>\n",
       "      <th>State</th>\n",
       "      <th>Zip</th>\n",
       "      <th>Bank</th>\n",
       "      <th>BankState</th>\n",
       "      <th>NAICS</th>\n",
       "      <th>ApprovalDate</th>\n",
       "      <th>ApprovalFY</th>\n",
       "      <th>...</th>\n",
       "      <th>RevLineCr</th>\n",
       "      <th>LowDoc</th>\n",
       "      <th>ChgOffDate</th>\n",
       "      <th>DisbursementDate</th>\n",
       "      <th>DisbursementGross</th>\n",
       "      <th>BalanceGross</th>\n",
       "      <th>MIS_Status</th>\n",
       "      <th>ChgOffPrinGr</th>\n",
       "      <th>GrAppv</th>\n",
       "      <th>SBA_Appv</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1000014003</td>\n",
       "      <td>ABC HOBBYCRAFT</td>\n",
       "      <td>EVANSVILLE</td>\n",
       "      <td>IN</td>\n",
       "      <td>47711</td>\n",
       "      <td>FIFTH THIRD BANK</td>\n",
       "      <td>OH</td>\n",
       "      <td>451120</td>\n",
       "      <td>28-Feb-97</td>\n",
       "      <td>1997</td>\n",
       "      <td>...</td>\n",
       "      <td>N</td>\n",
       "      <td>Y</td>\n",
       "      <td>NaN</td>\n",
       "      <td>28-Feb-99</td>\n",
       "      <td>$60,000.00</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>P I F</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>$60,000.00</td>\n",
       "      <td>$48,000.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1000024006</td>\n",
       "      <td>LANDMARK BAR &amp; GRILLE (THE)</td>\n",
       "      <td>NEW PARIS</td>\n",
       "      <td>IN</td>\n",
       "      <td>46526</td>\n",
       "      <td>1ST SOURCE BANK</td>\n",
       "      <td>IN</td>\n",
       "      <td>722410</td>\n",
       "      <td>28-Feb-97</td>\n",
       "      <td>1997</td>\n",
       "      <td>...</td>\n",
       "      <td>N</td>\n",
       "      <td>Y</td>\n",
       "      <td>NaN</td>\n",
       "      <td>31-May-97</td>\n",
       "      <td>$40,000.00</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>P I F</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>$40,000.00</td>\n",
       "      <td>$32,000.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1000034009</td>\n",
       "      <td>WHITLOCK DDS, TODD M.</td>\n",
       "      <td>BLOOMINGTON</td>\n",
       "      <td>IN</td>\n",
       "      <td>47401</td>\n",
       "      <td>GRANT COUNTY STATE BANK</td>\n",
       "      <td>IN</td>\n",
       "      <td>621210</td>\n",
       "      <td>28-Feb-97</td>\n",
       "      <td>1997</td>\n",
       "      <td>...</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "      <td>NaN</td>\n",
       "      <td>31-Dec-97</td>\n",
       "      <td>$287,000.00</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>P I F</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>$287,000.00</td>\n",
       "      <td>$215,250.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1000044001</td>\n",
       "      <td>BIG BUCKS PAWN &amp; JEWELRY, LLC</td>\n",
       "      <td>BROKEN ARROW</td>\n",
       "      <td>OK</td>\n",
       "      <td>74012</td>\n",
       "      <td>1ST NATL BK &amp; TR CO OF BROKEN</td>\n",
       "      <td>OK</td>\n",
       "      <td>0</td>\n",
       "      <td>28-Feb-97</td>\n",
       "      <td>1997</td>\n",
       "      <td>...</td>\n",
       "      <td>N</td>\n",
       "      <td>Y</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30-Jun-97</td>\n",
       "      <td>$35,000.00</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>P I F</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>$35,000.00</td>\n",
       "      <td>$28,000.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1000054004</td>\n",
       "      <td>ANASTASIA CONFECTIONS, INC.</td>\n",
       "      <td>ORLANDO</td>\n",
       "      <td>FL</td>\n",
       "      <td>32801</td>\n",
       "      <td>FLORIDA BUS. DEVEL CORP</td>\n",
       "      <td>FL</td>\n",
       "      <td>0</td>\n",
       "      <td>28-Feb-97</td>\n",
       "      <td>1997</td>\n",
       "      <td>...</td>\n",
       "      <td>N</td>\n",
       "      <td>N</td>\n",
       "      <td>NaN</td>\n",
       "      <td>14-May-97</td>\n",
       "      <td>$229,000.00</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>P I F</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>$229,000.00</td>\n",
       "      <td>$229,000.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 27 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   LoanNr_ChkDgt                           Name          City State    Zip  \\\n",
       "0     1000014003                 ABC HOBBYCRAFT    EVANSVILLE    IN  47711   \n",
       "1     1000024006    LANDMARK BAR & GRILLE (THE)     NEW PARIS    IN  46526   \n",
       "2     1000034009          WHITLOCK DDS, TODD M.   BLOOMINGTON    IN  47401   \n",
       "3     1000044001  BIG BUCKS PAWN & JEWELRY, LLC  BROKEN ARROW    OK  74012   \n",
       "4     1000054004    ANASTASIA CONFECTIONS, INC.       ORLANDO    FL  32801   \n",
       "\n",
       "                            Bank BankState   NAICS ApprovalDate ApprovalFY  \\\n",
       "0               FIFTH THIRD BANK        OH  451120    28-Feb-97       1997   \n",
       "1                1ST SOURCE BANK        IN  722410    28-Feb-97       1997   \n",
       "2        GRANT COUNTY STATE BANK        IN  621210    28-Feb-97       1997   \n",
       "3  1ST NATL BK & TR CO OF BROKEN        OK       0    28-Feb-97       1997   \n",
       "4        FLORIDA BUS. DEVEL CORP        FL       0    28-Feb-97       1997   \n",
       "\n",
       "   ...  RevLineCr  LowDoc  ChgOffDate  DisbursementDate  DisbursementGross  \\\n",
       "0  ...          N       Y         NaN         28-Feb-99        $60,000.00    \n",
       "1  ...          N       Y         NaN         31-May-97        $40,000.00    \n",
       "2  ...          N       N         NaN         31-Dec-97       $287,000.00    \n",
       "3  ...          N       Y         NaN         30-Jun-97        $35,000.00    \n",
       "4  ...          N       N         NaN         14-May-97       $229,000.00    \n",
       "\n",
       "   BalanceGross  MIS_Status ChgOffPrinGr        GrAppv      SBA_Appv  \n",
       "0        $0.00        P I F       $0.00    $60,000.00    $48,000.00   \n",
       "1        $0.00        P I F       $0.00    $40,000.00    $32,000.00   \n",
       "2        $0.00        P I F       $0.00   $287,000.00   $215,250.00   \n",
       "3        $0.00        P I F       $0.00    $35,000.00    $28,000.00   \n",
       "4        $0.00        P I F       $0.00   $229,000.00   $229,000.00   \n",
       "\n",
       "[5 rows x 27 columns]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv('SBAnational.csv')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "bad_columns = ['LoanNr_ChkDgt', 'NAICS', 'ApprovalDate', 'DisbursementDate', 'Name', 'FranchiseCode']\n",
    "target = 'MIS_Status'\n",
    "categorical = ['City', 'State', 'Zip', 'Bank', 'BankState', 'NewExist', 'UrbanRural', 'RevLineCr', 'LowDoc']\n",
    "ordinal = ['ApprovalFY', 'Term', 'NoEmp', 'CreateJob', 'RetainedJob']\n",
    "money_columns = ['DisbursementGross', 'BalanceGross', 'ChgOffPrinGr', 'GrAppv', 'SBA_Appv']\n",
    "date_col = ['ChgOffDate']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop bad columns\n",
    "df = df.dropna()\n",
    "df = df.drop(columns=bad_columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert date to ordinal days\n",
    "import datetime as dt\n",
    "for c in date_col:\n",
    "    df[c] = pd.to_datetime(df[c]).map(dt.datetime.toordinal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['CHGOFF' 'P I F']\n"
     ]
    }
   ],
   "source": [
    "# encode labels\n",
    "le = sklearn.preprocessing.LabelEncoder()\n",
    "df[target] = le.fit_transform(df[target].astype(str))\n",
    "class_names = le.classes_\n",
    "print(class_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# process ordinal features\n",
    "for col in ordinal:\n",
    "    df[col] = pd.to_numeric(df['ApprovalFY'].replace('1976A','1976', regex=False).replace('\\d+\\-\\w+\\-\\d+|,','', regex=True))\n",
    "    df[col] = df[col].astype('int32')\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# process money columns\n",
    "for c in money_columns:\n",
    "    df[c] = df[c].replace('\\$|,','', regex=True).replace('\\(','-', regex=True).replace('\\)','', regex=True)\n",
    "    df[c] = pd.to_numeric(df[c])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 1, 2, 3, 4, 8, 11, 12, 13]\n"
     ]
    }
   ],
   "source": [
    "# process categorical features\n",
    "cat_idxs = [df.columns.get_loc(c) for c in categorical if c in df]\n",
    "print(cat_idxs)\n",
    "categorical_names = {}\n",
    "for c in categorical:\n",
    "    le = preprocessing.LabelEncoder()\n",
    "    df[c] = le.fit_transform(df[c])\n",
    "    categorical_names[c] = le.classes_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = df[target]\n",
    "X = df.drop(columns=[target])\n",
    "\n",
    "# Split the data into train and test data:\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PDPs are defined as $\\hat{f}_{x_S}(x_S)=\\frac{1}{n}\\sum_{i=1}^n\\hat{f}(x_S,x^{(i)}_{C})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.10369109 0.05456078 0.11501452 0.05045366 0.03368958 0.01969201\n",
      " 0.01886796 0.021147   0.01493814 0.02170432 0.02113683 0.0109011\n",
      " 0.01481731 0.00260219 0.1321467  0.07635276 0.         0.17731801\n",
      " 0.04883301 0.06213304]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcYAAAEWCAYAAAD8XDcGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZhcRdn+8e9NWBIIEJFFCMIIoiwRAxNQVsOq8qoECYYYhIg/eFEUUUFwe2VRQVFRQIWgyBYg7AZEAhICCMEsZCNCBBKQTfYtIWzJ8/ujqslJp3umJ+mZ6Zm5P9fV13SfU6dO9WHIM3VO1VOKCMzMzCxZqbMbYGZm1kgcGM3MzAocGM3MzAocGM3MzAocGM3MzAocGM3MzAocGM3MzAocGM06gKRHJS2UNL/w2mgF6xws6Yl6tbHGc14o6Scdec5qJJ0k6dLObod1Pw6MZh3nsxHRt/B6qjMbI2nlzjz/iujKbbfG58Bo1skkfVzSPZJeljRD0uDCvi9LekDSa5LmSvrfvH0N4G/ARsUeaHmPrrxXmXuuJ0iaCSyQtHI+7hpJz0maJ+mYGtvdJClyGx+X9JKkoyTtIGlm/j7nFMqPlHS3pHMkvSLpQUl7FfZvJGmspBclPSzpiMK+kyRdLelSSa8CRwHfB4bl7z6jpetVvBaSviPpWUlPS/pyYX8fSb+S9Fhu3z8k9anhv9HIfK7X8vUbUcv1s8blv7rMOpGk/sBfgS8BNwN7AddI2jIingOeBT4DzAV2B/4maXJE3Cfp08ClEbFxob5aTjsc+B/geWAxcAPwl7x9Y+DvkuZExLgav8bHgC1y+8bm77E3sAowTdJVEXFHoezVwLrA54FrJX0gIl4ErgDuBzYCtgRulfRIRIzPx+4PHAQcCqyW6/hgRBxSaEvV65X3vw9YG+gP7ANcLen6iHgJ+CWwDbAz8N/c1sUt/TcCXgfOAnaIiDmSNgTWqfG6WYNyj9Gs41yfexwvS7o+bzsEuCkiboqIxRFxKzAF2A8gIv4aEY9EcgdwC7DbCrbjrIh4PCIWAjsA60XEKRHxVkTMBc4HDm5DfadGxBsRcQuwALg8Ip6NiCeBu4DtCmWfBX4TEW9HxBhgDvA/kt4P7AKckOuaDvyRFARLJkbE9fk6LazUkBqu19vAKfn8NwHzgQ9LWgk4HPhmRDwZEYsi4p6IeJNW/huR/rgYIKlPRDwdEbPbcO2sATkwmnWcIRHRL7+G5G2bAgcVAubLwK7AhgCSPi3p3nx78WXSP8brrmA7Hi+835R0O7Z4/u8DG7ShvmcK7xdW+Ny38PnJWHrlgsdIPcSNgBcj4rWyff2rtLuiGq7XCxHxTuHz67l96wK9gUcqVFv1v1FELACGkW7tPi3pr7knaV2YA6NZ53ocuKQQMPtFxBoRcbqk1YBrSLf4NoiIfsBNQOl+aaWlcRYAqxc+v69CmeJxjwPzys6/ZkTsV+G4euivpe/3bgI8lV/rSFqzbN+TVdq9zOcarldLngfeADavsK/qfyOAiBgXEfuQ/ph5kNTjti7MgdGsc10KfFbSJyX1ktQ7DxLZGFiV9CztOeCd/Exx38KxzwDvlbR2Ydt0YD9J60h6H3BsK+efBLyWB+T0yW0YIGmHun3Dpa0PHCNpFUkHAVuRblM+DtwDnJavwbbAV0jXp5pngKZ8GxRav15VRcRi4ALg13kQUC9JO+VgW/W/kaQNJO2vNBjqTdKt2cVtvCbWYBwYzTpRDgj7k25fPkfqnRwPrJRvKx4DXAm8BHyRNLildOyDwOXA3HyLbyPgEmAG8Cjp+dqYVs6/iDRYZSAwj9Rz+iNpgEp7+CdpoM7zwE+BoRHxQt43HGgi9R6vA34cEX9voa6r8s8XJN3X2vWqwXHALGAy8CLwc9J/h6r/jfLr27nNLwKfAL7ahnNaA5IXKjazjiBpJPD/ImLXzm6LWUvcYzQzMytwYDQzMyvwrVQzM7MC9xjNzMwKnBKuG1h33XWjqamps5thZtZlTJ069fmIWK/SPgfGbqCpqYkpU6Z0djPMzLoMSY9V2+dbqWZmZgUOjGZmZgUOjGZmZgUOjGZmZgUOjGZmZgUOjGZmZgUOjGZmZgUOjGZmZgWe4N8NTJ0KqmWNcjOzbqI903x3iR6jpPdJukLSI5KmSrpJ0pGSblyOuoZIminpAUmzJA0p7NtS0nRJ0yRtLumYXG60pJGSnsv7/yXpiCr1D5J0Vg3t2EDSZZLm5u80UdIBbf0+ZmZWXw0fGCWJtJr3hIjYPCKage8BGyxHXR8FfgnsHxFbAZ8Dfilp21xkCHB1RGwXEY8AXwP2iYgRef+YiBgIDAZ+JmmDsvpXjogpEXFMDd/peuDOiNgsf6eDgY0rlHWv3sysA3WFf3T3AN6OiHNLGyJihqT3AHtJuhoYAEwFDomIkLQf8GtgAXA3sFlEfAY4DvhZRMzL9cyTdBpwvKTLgWOBRZL2AuYAmwF/k3QB8FLh/M9KegTYVNJXgc1z2f9IOg84LiI+I+kkYJO8bxPgNxFxFrAn8FbZd3oMOBveXen880BfoBfwiXpdTDMza1lXCIyloFfJdsA2wFOkALiLpCnAecDuOfBdXii/DanHWDQFODoibpJ0LjA/In4JIOlTwB4R8XwOVuTtm5GC3cN509bArhGxUNLgsvq3JAX3NYE5kv6Q23FfK997e2DbiHix0k5JRwJHpk+btFKVmZnVquFvpbZiUkQ8ERGLgelAEykQzS31CoHLqx28HIZJmp7r/N9C0BobEQurHPPXiHgzIp4HnqXCLWBJv5M0Q9LkwuZbqwVFgIgYFRGDImIQVFw5xczMlkNXCIyzgeYq+94svF9E6z3gf1WoqzmfoxZjImJgRHwsIq4rbF/QwjGV2jib1CMEICKOBvZi6QjXUp1mZtZOukJgHA+slm8dApAHy+xWpfwcYDNJTfnzsMK+XwLfK+3LP78P/KqO7a3FeKB3fj5ZsnoHt8HMzCpo+GeMeTDNAcBvJJ0AvAE8ShrVWan8QklfA26WtACYXNg3Pddxg6RVgLeB70bE9Pb+HmVtjDxN5ExJ3wWeI/UQT1ie+pqbwesUm5nVh6I9Z0l2Ekl9I2J+nhbxO+ChiDizs9vVXgYNGhRTHBnNzGomaWoao7Gshu8xLqcjJB0GrApMI41S7bac+casa+mG/ZFupSs8Y2xVeWYcYB/g96TRqSMi4vU21LWimXGmSXpI0jhJO9d4vq2X64ubmVnddfnA2ICZcbaLiC2A04FrJW3VymmHkOZBmplZA+jygZEqmXGAu4C+kq6W9GDu1QlA0n5521RJZxVyri6TGQcoZcbZj5QZ56uSbs/JAEqZcb5V3qiIuB0YRZ6EL+kISZPzfMVrJK2ee5SfA87IPdHN8+vm3La7JG3ZLlfNzMwq6g6BsbXMOMeSemSbkTLj9CY9c/x07l0W5w5uU6GuKcA2EXETcC5wZkTsERFHkTLu7NHCwJ77SAkHAK6NiB0i4qPAA8BXIuIeYCxwfJ4f+QgpmH4jt+040i3hZSglUZ+SMv08V+X0ZmbWVt118E3JpIh4AiBnrGkC5rNsZpwjKx++wopDYgZI+gnQj5QDddwyhaW+wM7AVVoymma1ShVHxChSEEUa5Ef5ZmZ10h0C42xgaJV9y5sZZ0ZhW1sy45TbjtQ7BLgQGJIToI8krdBRbiXg5byCh5mZdYLucCu1ITPjSPoEqSd6ft60JvB0TiwwolD0tbyPiHgVmCfpoFyH8oAgMzPrIF2+x9hgmXGGSdqVlN5tHnBgRJR6jD8C/kl6IPhPcjAErgDOl3QMqec7AviDpB8Cq+T9xR7sMpz5xsysfrpl5pvWdLfMOM58Y2bWNj0x801rulVmHGe+MWu7HtgnsBr1yB5jZ8i3e39ctnlb4EvAARFRbQBRDXUPijSrxMxq5X/6eraWeowOjJ0kDxYaQZoHuXjF6nJgNGsr/9PXs7UUGLvDqNQuR9KHgP8j9RY3kXR/3j5S0l8kTcj5Vst7mGZm1s4cGDtYHu16GfCdiPhPhSI7AgeSbrMeJKlyV9+Zb8zM2oUDY8c7FZgdEWOq7L81Il6IiIXAtcCulQpFxKiIGJRuBaxXqYiZmS2HnjoqtVNIGkzqDW7fQrHyJx9+EmJm1oHcY+wgkt4D/Bk4NCJea6HoPpLWkdSHtCTV3R3SQDMzA9xj7EhHAeuTstoUt19eVm4ScA2wMXBpRLQ63NSZb8zM6seBsYNExGmktR0r+Xnh/RMRMaQDmmRmZhU4MHYDznxjPYHnHVpH8TPGdibpfZKukPSIpKmSbpK0u6Sr8/6BkvYDiIgLI+LrndtiM7OezYGxHeUk5dcBEyJi84hoBr5HWhSklAJuILBfZ7XRzMyW5sDYvvYA3o6Ic0sbImIG8Lik+yWtCpxCWq5quqRhOePNegCSVpL0cOmzmZm1PwfG9jUAmFptZ0S8RUoNNyYiBuZJ/5eyZCHjvYEZEbFMahtnvjEzax8OjI3nAuDQ/P5w0tzHZTjzjZlZ+3BgbF+zgea2HBARjwPPSNqTlDf1b+3RMDMzq8yBsX2NB1bLS0wBIGlb4P2FMq8Ba5Yd90fSLdWrImJRu7fSzMze5cDYjiItdnkAsHeerjGbNMn/v4VitwNblwbf5G1jgb5UuY1arrk5zfHyy6/u/DLrKJ7g384i4ingCxV2Dcj7XwR2KNv3UdKgmwfbuXlmZlbGgbHBSDoR+CpLRqa2yplvrDXucZnVrtVbqZIW5dt8syXNkPQdSSvlfYMkndXCsYMl3VjPBncWSd8v+7yBpMskzc0ZbSZKOmBFzxMRp0fEphHxjxWty8zM2q6WZ4wL8xy7bYB9gE8DPwaIiCkRcUx7NU5SI/Vo3w2MOaPN9cCdEbFZzmhzMGlFjKU02HcwM7NWtGnwTUQ8CxwJfF3Juz1CSZ/IPcvpkqZJKo20XEvSXyXNkXRuobc5v1SvpKGSLszvL8zl/gn8olq9ko6XNFnSTEkn521Nkh7Mdfxb0mhJe0u6O2eU2TGXW0PSBZIm5Tr3z9tHSrpW0s25/C/y9tOBPrkNo4E9gbfKMto8FhFnF+oZK2k8cFu+VmfkbDezSoNsJG0o6c5c7/2SdpPUK7e/VPZbbflvZGZmK6bNvZmImCupF2ltwaLjgKMj4m5JfYE38vYdga2Bx4Cbgc8DV7dymo2BnSNikaQbyuuVtC+wRa5bwFhJuwP/AT4IHESaHD8Z+CKwK/A5Uq9vCPADYHxEHC6pHzBJ0t/zuQcC2wFvAnMknR0RJ0r6ekQMBJB0DHBfK99he2DbiHhR0oG53o8C6wKTJd2Z2zYuIn6ar+nquVz/iBiQz9WvUuV5CkieBrJJK00xM7Na1XO6xt3Ar3PQ6BcR7+TtkyJibp6PdzkpSLWmOH+vUr375tc0UoDakhQoAeZFxKyIWEyaYH9bnjYxC2jKZfYFTpQ0HZgA9GZJdLktIl6JiDeAfwGbttZYSb/Lz18nFzbfmkeckr/z5RGxKCKeAe4gjUSdDHxZ0knARyLiNWAusJmksyV9Cni10jmd+cbMrH20OTBK2gxYBDxb3B4RpwP/D+gD3C1py9KusiqiwvbeZWUWtFKvgNPys8+BEfHBiPhTPuTNQj2LC58Xs6SHLODAwvGbRMQDFY5fROVe9WxSj7DUxqOBvVg6Qi0oP6hcRNwJ7A48CVwo6dCIeInUs5wAHEWa7G9mZh2kTYFRaZWHc4Fzci+suG/z3FP7OaknVAqMO0r6QH62OAwojbZ8RtJWeXvV0ZxV6h0HHJ5vrSKpv6TyW7stGQd8Iw+iQdJ2NRzztqRV8vvxQG9JXy3sX72FY+8iraDRK1/D3Um3bzcFnomI80kBcHtJ6wIrRcQ1wA8pBGAzM2t/tTxj7JNvOa4CvANcAvy6QrljJe1B6pnNJuX43IkUzM4hPfu7nbQ+IcCJwI2kpSGmkDK9VLJMvRHxpqStgIk5ts0HDiH18GpxKvAbYGYOzPOAz7RyzKhc/r6IGCFpCHCmpO/m77AAOKHKsdeRrsUMUk/5uxHxX0mHAcdLejt/h0OB/sCfS4OUSOs3tqi5GaZMaa2UmZnVQuGZv13eoEGDYoojo5lZzSRNTWM0luU5dt2AM990Xf671KzxOIm4mZlZQY8IjJJ+oJTSbmaeTP8xSROUkg5Ml/SACktD5WMGSoo8ZaKWcwzJ5bdsvbSZmTWqbh8YJe1EGlizfURsC+wNPJ53j8iT9ncBfi5p1cKhw0kjaIfXeKq2ljczswbU7QMjsCHwfES8CRARz+eloIr6kkaVLoJ3c6EeBIwE9pFUPs9yKXnayK7AV0g5U0vbB+eUbxVT4kk6M/dkb5O0nqQtJU0qHN8kaVaVcx4paYqkKWlQrJmZ1UNPCIy3AO9Xyp36e0mfKOwbLWkmMAc4tZBtZ2dSBp1HSBPt/6eVc+wP3BwR/wZekNRc2Lcj8A1SWrzNSSnxANYApuTk7HcAP87rL64q6QO5zDBgTKUTOvONmVn76PaBMSLmA82kvKLPAWMkjcy7R+Tbq5sAx+UJ95Buh16R319B67dHWypfLSXeYpYEvUsL268kBURoITCamVn76BHTNXJQmgBMyLcmDyvb/5yk+4CPSXoCOBDYX9IPSOnj3itpzZzLdCmS1iGttvERSQH0AkLS8aXqy5tTrZn55xjgKknXpqbFQ238umZmtgK6fY9R0oclbVHYNJC00kexzOqkFTUeIeU8nRkR74+IpojYFLiG6mnrhgKX5MWFmyLi/aRMOrvl/dVS4q2Uj4W0ysY/APLt20XAj6ixt9jcnObD+dX1XmbWeLp9YCQNrLlI0r/y88StgZPyvtE53d1U4MKImEq6DXpdWR3XUP12amvlSynxHiAFzFLZBaSgeT+px3lK4fgxpBR3V9b4Hc3MrE6cEq4dSRoMHBcRy+RhlTQ/Iqrlh23jeQZFSjdrncn/K5l1HS2lhOsJPcZ2VyWBwLHAajUce2y+lWtmZg2gRwy+qQdJ7wVuq7DrBJYkEHgzLxu1Kul26KBKvUWAQm/xWNKo1Nfr32ozM2sr9xhrFBEvFBY2fvdFmo+4VAIB0qCajYDbJd0OIOkPeUL+bEkn523HVCi3r6SJku6TdFVpzUkzM+sYfsa4gnLg+gdpoeK/A2Mi4g5Jj5J6jM/ncutExIuSepF6nsdExMxiudzbvBb4dEQskHQCsFpEnFLhvEeS5mYCmzSXDbS1TuD/lcy6Di871Y4iYn7OdLMbsAcpgcCJFYp+IQezlUlp6rYGZpaV+XjefndegHlVYGKV844iLZ6cB9+YmVk9ODDWQWsJBHKKt+OAHSLiJUkXApXyrwq4NSKciNzMrJP4GeMKaiGBwGvAmnnbWqR5i69I2gD4dKF8sdy9wC6SPpjrXkPSh9qz/WZmtjT3GFdcX+BsSf2Ad4CHSc/+hgM3S3oqIvaQNA14kLTk1d2F40eVlRsJXC6pNNXjh8C/W2pAczNM8TRGM7O68OCbbmDQoEExxZHRzKxmHnzTzU2dCmmsjrUn/w1p1jP4GWMdSFqUM97MyPMPd16BuiZIqvhXjJmZtT/3GOtjYZ7sj6RPAqcBn+jcJpmZ2fJwj7H+1gJegjT5X9JtuRc5S9L+eXuTpAcknZ8z4dwiqU+xEkkrSbpQ0k864TuYmfVY7jHWR5+8fFVv0uT9PfP2N4ADIuLVnNXmXklj874tgOERcYSkK0mLI1+a960MjAbuj4ifVjphWeab+n8jM7Meyj3G+liYc6duCXwKuFgpdY2An+V1IP8O9Ac2yMfMi4jp+f1UoKlQ33m0EBQhZb6JiEFpVNV6df46ZmY9lwNjnUXERGBdUrQakX8252eQz7Ak482bhcMWsXTv/R5gD0mVsuOYmVk7cmCsM0lbAr2AF4C1gWcj4m1JewCb1ljNn4CbgCsl+Xa3mVkH8j+69VF6xgjp9ulhEbFI0mjghpw/dQop801NIuLXktYGLpE0IiIWVyvrzDdmZvXjzDfdgDPfmJm1jTPfdHPOfNN+/HejWc/jZ4xVSNpA0mWS5kqaKmmipANaKH+spDfy7U8zM+uiHBgryFMtrgfujIjNIqIZOBjYuKxcscc9HJgMfL7DGmpmZnXnwFjZnsBbEXFuaUNEPBYRZ0saKWmspPHAbQCSNictP/VDUoAkbx8p6S85/+lDkn6ctzdJelDS6JwB52pJq0v6lKSrCscPlnRjB31nMzPDgbGabYD7Wti/PTA0Ikr5UA8GrgDuAj6cFyMu2ZGU1WZb4KBCgvAPA7+PiK2AV4GvkZIAfEzSGrnMsFzvMiQdKWmKpCnwXJu/oJmZVebAWANJv8srZ0zOm26NiBcLRYYDV+QpFdcABxX23RoRL0TEQuBaYNe8/fGIKC1YfCmwa0S8A9wMfDbfpv0f4C+V2uTMN2Zm7cOjUiubTerlARARR+dcp6U5EQtK+yR9hJT39Nb0aJJVgXnAOaXDy+qOVrZfAXwdeBGYEhGvrdA3MTOzNnGPsbLxQG9JXy1sW71K2eHASRHRlF8bARtJKmW52UfSOnn1jCFAqZe4iaSd8vsvAv/I7+8g3ao9giq3Uc3MrP04MFYQKevBEOATkuZJmgRcBJxQofjBwHVl267L2wEmkW6vzgSuiYhSr3MOcLSkB4D3AH/I514E3Ah8Ov9sVXNzmm/nV/1fZtbzOPNNO5I0EhgUEV8v294E3BgRA+pxHme+MTNrG2e+6eZ6WuYb/y1nZu2pW99KlbRI0vQ8ovQ+STuvQF0TClMtitsPlzRL0kxJ90vaP28fCdxS3lsEiIhHS73FPNdxo+Vtl5mZ1Vd37zEuzOsgIumTwGnAJ1o+pHaSNgZ+AGwfEa9I6suSuRMjgfuBp1qpptZyZmbWAbp1j7HMWsBLAJL6Srot9yJnFXp5TTkTzfmSZku6JY8mfZeklSRdKOknwPrAa8B8gIiYHxHzJA0FBgGjc4+1j6T/kzQ59ypHKalUrlnSHTk/6zhJG3bcJTIzs+4eGPvkgPMg8Efg1Lz9DeCAiNge2AP4Vc6PCmlO4u8iYhvgZQrzGUk97NHAQxHxQ2AG8AwwT9KfJX0WICKuJs15HBERA/Pk/nMiYod8C7UP8JnycsA7wNmkrDrNwAXATyt9MWe+MTNrHz3pVupOwMWSBpAWE/6ZpN2BxUB/oJTGbV5ElBYdngo0Feo7D7gyIn4KaWqFpE8BOwB7AWdKao6Ikyq0ZQ9J3yXNh1yHlETghrIyHwYGsCRZQC/g6UpfLCJGAaPSdxvk4ShmZnXS3QPjuyJiYs5esx6wX/7ZHBFvS3oU6J2Lvlk4bBGpd1dyDynA/Soi3sj1Bmmu4iRJtwJ/Bk4qnltSb+D3pKkbj0s6qXC+pYoCsyNipwr7zMysA3T3W6nvkrQlqQf2ArA28GwOinsAm7Z48BJ/Am4CrpS0sqSNJG1f2D8QeCy/fw1YM78vBcHn8wCdoYVjiuXmAOuVMuJIWkXSNjV/STMzW2HdvcfYR1LptqiAw/Ltz9HADZJmkZ7xPVhrhRHxa6XFiC8BTgR+madbvEF62HdULnohcK6khcBOwPmk0af/Ja3bSJVyQ4Gz8jlWBn5Duu1aVXMzeH6/mVl9OPNNN+DMN2ZmbePMN91cI2e+8d9dZtbV9JhnjACS3ifpCkmP5HmCN0n6UB3qPVZStdU3iuUqZs8pKzN/RdtjZmbLr8cExjxP8TpgQkRsnucJfo8l0zTIiwMvj2OpviyVmZl1IT0mMJIm8r8dEeeWNkTEDKCXpLskjQX+BSDpEEmTcnKA8yT1ytv/kCfVz5Z0ct52DLARcLuk2/O2fSVNzJl1rsojUZciaXjOunO/pJ+X7Tszn+M2SeuVH2tmZu2nJwXGAaQJ+5VsD3wzIj4kaStgGLBLTg6wCBiRy/0gP6zdlrRW47YRcRYpz+keEbFHniv5Q2DvnFlnCvDt4snyKNafA3uSpnjsIGlI3r0GMCVn3rkD+HGlBjvzjZlZ+/Dgm2RSRMzL7/cCmoHJOftMH+DZvO8Lko4kXbcNga1JCxAXfTxvvzsfvyowsazMDqRbus8B5OkjuwPXkzLxjMnlLgWurdRgZ74xM2sfPSkwzmbpifVFCwrvBVwUEd8rFpD0AeA4YIeIeEnShVTPXnNrRAxf8SYD4KBnZtaBetKt1PHAarnHB4CkbYHdysrdBgyVtH4us46kTUmrcywAXpG0AfDpwjHF7DX3ArtI+mA+fo0KI18nkW7FrpufXw4n3TaF9N+kFMC/CPxjeb+wmZm1XY8JjDmn6QHA3nm6xmzS+oz/LSv3L9IzwlskzQRuBTbMA3WmkbLkXAbcXThsFHCzpNvz7dGRwOX5+InAlrncysCbEfE0KWvO7aQVOqZGxF9ymQXAjpLuJz2DPKW179bcnOYLNuLLzKyrceabDiJpNeBhYEBEvFLPup35xsysbZz5ppPlSf2XAL+vd1CExs1847+5zKwrcmBsgaQAfh0R38mfjwP6VllvsXTMScARLD2HYnBEbNWOTTUzszrpMc8Yl9ObwOfz3MS2ODMiBhZeL7dH48zMrP4cGFv2DmlgzbfKd0hqkjRe0sycoWaTliqSNFLS9ZJulfSopK9L+rakaZLulbROLjdB0m9z1p37Je3YPl/NzMwqcWBs3e+AEXl9xKKzSfMdtwVGA2cV9n0rB7bppTRx2QDg86QJ/j8FXo+I7UgjVw8tlFs9Z935GnBBpUY5842ZWftwYGxFRLwKXAwcU7ZrJ9K0DUgDa3Yt7CveSt2jsP32iHgtT+l4Bbghb58FNBXKXZ7PfSewlqR+Fdo1KiIGpVFVTqdqZlYvDoy1+Q3wFVIe0xXxZuH94sLnxSw9EKp8PKfHd5qZdRAHxhpExIvAlaTgWHIPcHB+PwK4q46nHAYgaVfglfaY4mFmZpV5ukbtfgV8vfD5G8CfJR1Pesj35cK+b0k6pPB5CG3zhqRpwCrA4a0Vbm4Gz+83M6sPZ75pMJImAMdFRM2hzplvzMzaxplvujlnvjEzq58u/YxR0qLCfL8bKo3eLCvfT9LXaqz7njq1sSknBEfSYEk3tnLIBGBwPc5tZmZt16UDI7AwT4kYALwIHN1K+X6kuYGtioidVzA7Ki0AABtXSURBVLRxZmbW9XT1wFg0Eehf+iDpeEmTc2aak/Pm04HNcy/zDEl9c9aa+yTNkrR/4fj5+efgnI3makkPShotpRuXkpol3SFpqqRxkjYsbJ8haQZVgnVe5/H63L5789qQJR+VNFHSQ5KOqOtVMjOzFnWLZ4x5sd+9gD/lz/sCWwA7AgLGStqdtAbigJxVBkkrAwdExKs5H+q9ksbGsiOStgO2AZ4ircO4i6R/krLf7B8Rz0kaRspmczjwZ+DrEXGnpDOqNPtkYFpEDJG0JymJwMC8b1vg46R5k9Mk/TUinir7zkcCedHlFrPRmZlZG3T1wNhH0nRST/EB0qLCAPvm17T8uS8pUP6n7HgBP8tBc3GuZwPKFi8GJkXEEwD5fE3Ay6QUb7fmDmQv4On8nLNfzloDKSvOpyu0fVfgQICIGC/pvZLWyvv+EhELgYU5pdyOwPXFgyNiFCmPK9IgD3MxM6uTrh4YF0bEQEmrA+NIty3PIgW80yLivGJhSU1lx48g5VNrjoi3JT0K9K5wnmLGmkWk6yZgdkTsVHaOFgcA1ciZb8zMOkm3eMYYEa+Tcpl+J98eHQccLqkvgKT+ktYHXgPWLBy6NvBsDop7AJu24bRzgPUk7ZTPsYqkbfISUy/nrDWQgm8ld5X2SRoMPJ/zsgLsL6m3pPeSRqhObkO7zMxsBXT1HuO7ImKapJnA8Ii4RNJWwMR8m3M+cEhEPCLp7jx94m/Az4EbJM0CpgAPtuF8b0kaCpyVV95YmZRTdTYpC84FSgsd31I4bGWW9D5PymVmAq8DhxXKzQRuB9YFTi1/vljOmW/MzOrHmW86kKRvAv0j4rv1rNeZb8zM2saZbxqApD+RBut8od51N1LmG/+dZWZdXbd4xlgkaYikkLRlZ7elKCK+Qrp9exCApAslzSssaHxMniP51dIxkj6W5zmu0lntNjPrabpjj3E48I/888crWpmklSPinRVuVWXHR8TVhXNtQHouejXwAnAO8LWIeLudzm9mZmW6VY8xj0LdlbRu4sF522BJd0r6q6Q5ks6VtFLeN1/SmZJm5ww46+XtEyT9RtIU4Js53+n43Hu7TdImktaW9FihrjUkPZ5Hpx6Rs+7MkHRNnk7Sqoh4Bvgl8AvgKGBmRPyj3tfJzMyq61aBEdgfuDki/g28IKk5b9+RtH7i1sDmwOfz9jWAKRGxDXAHS/cwV42IQRHxK1KGm4siYltgNHBWXjx4OvCJXP4zwLjcu7s2InaIiI+SEg8UFzguOqNwK/Ujedu5uZ3HA1UH6Ug6UtKUFLyfq+XamJlZDbpbYBwOXJHfX5E/Q8pcMzciFgGXk3qVkLLdjMnvLy1sp7AdYCfgsvz+kkK5McCw/P7gwjEDJN2Vp4GMIKWTq+T4nAR9YETMAoiIxcB5wN8i4oVqXzQiRuXAPSjlKDAzs3roNs8YJa0D7Al8JM8f7EXKGPNXas8kU9y+oIbTjiWllFsHaAbG5+0XAkMiYoakkbR9GanF+WVmZh2sO/UYhwKXRMSmEdEUEe8H5gG7ATtK+kB+HjiMNDgH0vcfmt9/sbC93D3kZ5akHuBdABExn5SV5rfAjblHCim7ztN5NGm1zDdmZtaAulNgHA5cV7btmrx9MmmE5wOkYFkqt4AUNO8n9TZPqVL3N4Av5yw1XwK+Wdg3BjiEpW+9/gj4J2kljpqz6Syv5uY0f7ARXmZmXV23z3yT85AeFxGfqbBvfkT07fhW1Zcz35iZtY0z33RzjZL5ppv/jWVmPUS3D4wRMQGYUGXfCvcW8woYt+WP7yMtS1WaP7FjRLy1oucwM7OO0+0DY3vLUyoGAkg6CZgfEb+s9XhJvQqDdszMrJN1p8E3DUfSYZIm5Qn8v5e0kqSVJb2cM+vMJA3+eULSz3KmnMmStpd0i6RHJB3R2d/DzKwncWBsJ5IGAAcAO0fEQFLvvDTlY23gzojYNiIm5m3zcqace4E/lY4FTq1SvzPfmJm1A99KbT97AzsAU/JiyX2Ax/O+t1h2asnY/HMWsHJELAAWSFosqW+eM/muiBgFjAKQBnnYi5lZnTgwth8BF0TEj5baKK0MLIxl58m8mX8uLrwvffZ/JzOzDuJbqe3n78AXJK0LafSqpE06uU1mZtYKB8Z2kpOCnwz8PQ+yuQXYoD3O1SiZb8zMuoNun/mmJ3DmGzOztnHmm26uozLf+G8oM+sJeuStVEkh6VeFz8flyfkrUudgSa8UFh6eLmnvFspvJOnqFvb3k/S1FWmTmZm1XY8MjKRRn58vDYypo7sKCw8PjIi/VysYEU9FxNBq+4F+gAOjmVkH66mB8R3SHMBvle+QtJ6ka3IGmsmSdsnbZ+VenCS9IOnQvP1iSftUO5GkHSTNlNRb0hqSZksaIKkpL3eFpG0KGXJmStoCOB3YPG87oz0ugpmZLasnP2P8HTBT0i/Ktv8WODMi/pGnV4wDtiKtrbgL8Bgwl7QA8sXATsBXSZP5d5M0vVDXgRExWdJY4CekSf6XRsT9kpoK5Y4CfhsRoyWtCvQCTgQG5Kw5y5B0JHBk+uRZIGZm9dJjA2NEvCrpYuAYYGFh197A1loymmUtSX2Bu4DdSYHxD8CRkvoDL0XEglz+rkrrPpIWQJ4MvJHPV24i8ANJGwPXRsRDamU0jTPfmJm1j556K7XkN8BXgDUK21YCPl54Ttg/p2O7k9RL3I20jNVzwFBSwGzNe4G+wJpA7/KdEXEZ8DlSgL5J0p7L/Y3MzGyF9OjAGBEvAleSgmPJLcA3Sh8kDcxlHwfWBbaIiLnAP4DjSAGzNecBPwJGAz8v3ylpM2BuRJwF/AXYFniNFEjNzKwD9ejAmP2KFPBKjgEG5UEw/yI9/yv5J/Dv/P4uoD8pQJbsVjZdY2gepPN27hWeDuxQoUf4BeD+/HxyAHBxXufxbkn3tzb4pqMy35iZ9QTOfNMNOPONmVnbOPNNN+fMN2Zm9dPtb6VKWpRva94v6QZJ/ZaznnfnHZZtP6WlDDc11LujpDslzZE0TdIfJa2+vPWZmdmK6faBkbT24cCIGAC8CBxdz8oj4v9aynDTEkkbAFcBJ0TEhyNiO+Bmygbd5DUczcysA/SEwFg0kTRgBgBJx+fsNjMlnZy3nS7p6EKZkyQdV61CSRdKGprfPyrpZEn35Uw5W+bta0i6IGe3mSZp/3z40cBFETGxVF9EXB0Rz+TzXiLpbuCSel4EMzOrrscERkm9gL2AsfnzvsAWwI7AQKBZ0u7AGNIo0ZIv5G21ej4iticlASgF1B8A4yNiR2AP4AxJa5BGoE5toa6tgb0jYniF73OkpCmSpqQplWZmVg89ITD2ydMg/ktaKPjWvH3f/JoG3AdsSZqjOA1YP69+8VFSZpvH23C+a/PPqUBT4Vwn5nZMIE3yryWP29iIWFhpR0SMiohBaVTVem1onpmZtaQnPLtaGBED84CWcaTbl2cBAk6LiPMqHHMVKavN+2hbbxHSyh0Ai1hyfUXKmzqnWFDSbKCZNKm/kgVtPLeZma2gntBjBCAiXidN3v9OHswyDjg850FFUn9J6+fiY4CDScHxqjqcfhzwDeUEqJK2y9vPAQ6T9LFSQUmfz4NyzMysE/SYwAiQb5POBIZHxC3AZcBESbOAq8mjQSNidn7/ZEQ8Xajiw5KeKLwOqvHUpwKrkFbzmJ0/ExHPkALwL/N0jQeAT5LSwdXMmW/MzOrHmW+6AWe+MTNrG2e+6ebqmfnGfyeZWU/XbW6lVspMU20OYnHuYR3OOyHfBp2R50RWXFh4OeseLOnGetVnZmat6zaBsVbtlEVmRER8FPg90OJKGB3UHjMzW049IjDmXt1v0mR4vpk3750nyP9b0mdyuSZJd+XMNfdJ2jlvH5zruFrSg5JGl0aYlinPrDO/8H6opAvz+wslnSvpn8Avcr7UiTkrzj2SPtw+V8LMzFrTk3orq5YetOYA1UTKerM5cLukDwLPAvtExBuStgAuB0oPZ7cDtgGeAu4GdmHptRgBPgVcX2N7NgZ2johFktYCdouId3JC8p8BB7Z0sKQjgSPTp1pyBZiZWS26U2CsNmyktL18ov6VEbEYeEjSXFLmm3nAOfk54SLgQ4XykyLiCYCcwaaJJYFxtKRVgb6k9HK1uCoiFuX3awMX5WAcpKkdLYqIUcCo1J5BHjJjZlYn3elW6gvAe8q2rQM8n9+XZ5EpDyYBfAt4Bvgoqae4amH/m4X3xaw2ACOAzYCLgLOrnKN32fmK7TkVuD2vAPLZCmXNzKyDdJvAGBHzgacl7QkgaR3Src3y250lB0laSdLmpKA2h9Rzezr3JL8E9GrD+QP4EfDx0qoawDOStpK0EnBAC4evDTyZ34+s9ZxmZlZ/3SYwZocCP8q3OscDJ0fEI1XK/geYBPwNOCoi3iCNKj1M0gzSrdU25SrNCb9/BRyfN50I3AjcAzxd7TjgF8BpkqaxHLe365n5xsysp3Pmm27AmW/MzNrGmW+6ufLMN/5bx8xs+XW3W6krrDj3sA51jZT0XJ6f+JCkcaW5kWZm1pgcGNvfmIjYLiK2AE4HrpW0VWc3yszMKnNgrEHOiDNe0kxJt0naRFIvSfOU9JO0SNLuufydeU7iUiLidtLcwyNzuYGS7s31XifpPXn7ByX9PedfvS+PnDUzsw7gwFibs4GLImJbYDRwVp6cPwfYGtgVuA/YTdJqwPsj4qEqdd1HGvEKcDFwQq53FvDjvH008Lucf3VnKoxolXRkTmk3BZ6ry5c0MzMHxlrtRFrUGOASUiAEuAvYPb9Oy9t3ACa3UJcAJK0N9IuIO/L2i4DdJa0J9I+I6wAi4o2IeL28kogYFRGD0qiq9Vboy5mZ2RIOjCvmTmA3Us7Vm4B+wGBSwKxmO+CBdm+ZmZktFwfG2twDHJzfj2BJ4JtEutW5OCcImA78LylgLkPSJ0jPF8+PiFeAlyTtlnd/CbgjIl4DnpA0JB+zmqTV2+E7mZlZBQ6My1pd0hOF17eBbwBfljSTFMC+CRARbwKPA/fmY+8C1iQ9LywZJmm6pH8D3wcOjIhSj/Ew4Ixc70DglLz9S8Axefs9wPtaanB55hszM1t+znzTDTjzjZlZ27SU+cY9RjMzs4KGDYx5XuD0wly+VjPG1DNrTQ3n2kDSZZLmSpoqaaKkllbQMDOzLqBhAyOwMCIG5rl83yNNh2gIkgRcD9wZEZtFRDNpcM7GFco6H62ZWRfSyIGxaC3gJQBJfXP2mfskzZK0f3nhamVyBpsHJJ0vabakWyT1yfsqZpuRdLykyTk7zcn5FHsCb0XEuaVzRsRjEXF2PmakpLGSxgO35ew4Z0i6P7dnWC63Yc6SMz3v2y1n1LmwUPZb7XdZzcysXCP3ZvrkdRV7AxuSghHAG8ABEfGqpHWBeyWNjaVHEVUsk/dtAQyPiCMkXQkcCFxKyjZzekRcJ6k3sJKkfXP5HUkT88fmtG/bkDLYtGR7YNuIeFHSgaRRpx8F1gUmS7oT+CIwLiJ+KqkXsHou1z8iBgBI6lepcklHklPLbbLJJq00xczMatXIgXFhRAwEkLQTcLGkAaQA9bMcoBYD/YENgP8Wjq1WBmBeREzP76cCTZWyzeTz7gvsC0zL5fuSAuVSJP2OlPXmrYjYIW++NSJezO93BS7PaeSekXQHSzLkXCBpFeD6iJguaS6wmaSzgb8Ct1S6OBExipR3lUGDBnlosZlZnXSJW6kRMZHU01qPNMF+PaA5B85nSL3KopbKvFkot4iW/zgQcFp+1jkwIj4YEX8CZpN6hKX2HQ3sxdK52RbU8L3uJKWTexK4UNKhEfESqWc5ATgK+GNr9ZiZWf10icAoaUugF/ACsDbwbES8LWkPYNMKh9RS5l0tZJsZBxwuqW/e3l/S+sB4oLekrxaqaSk7zV2kif69JK1HCoaTJG0KPBMR55MC4Pb51u9KEXEN8EMKAdjMzNpfI99KLT1jhNRzOywiFkkaDdwgaRYwBXiwwrG1lCn3JeA8SacAbwMHRcQtSmsnTkwDUZkPHBIRz+Ygeqak75KWt1gAnFCl7utIichnAAF8NyL+K+kw4HhJb+e6DyXd9v2zpNIfLd+roe1mZlYnznzTDTjzjZlZ2zjzjZmZWY0cGM3MzAocGM3MzAocGM3MzAocGM3MzAocGM3MzAocGM3MzAocGM3MzAo8wb8bkPQaMKez29FG6wLPd3Yj2qCrtRfc5o7iNneMerd504hYr9KORk4JZ7WbUy2DQ6OSNKUrtbmrtRfc5o7iNneMjmyzb6WamZkVODCamZkVODB2D6M6uwHLoau1uau1F9zmjuI2d4wOa7MH35iZmRW4x2hmZlbgwGhmZlbgwNhgJH1K0hxJD0s6scL+1SSNyfv/KampsO97efscSZ+stc7OarOkfSRNlTQr/9yzcMyEXOf0/Fq/QdrcJGlhoV3nFo5pzt/lYUlnSVKDtHlEob3TJS2WNDDv6+zrvLuk+yS9I2lo2b7DJD2UX4cVtnf2da7YZkkDJU2UNFvSTEnDCvsulDSvcJ0HNkKb875FhXaNLWz/QP49ejj/Xq3aCG2WtEfZ7/MbkobkffW5zhHhV4O8gF7AI8BmwKrADGDrsjJfA87N7w8GxuT3W+fyqwEfyPX0qqXOTmzzdsBG+f0A4MnCMROAQQ14nZuA+6vUOwn4OCDgb8CnG6HNZWU+AjzSQNe5CdgWuBgYWti+DjA3/3xPfv+eBrnO1dr8IWCL/H4j4GmgX/58YbFso1znvG9+lXqvBA7O788FvtoobS77PXkRWL2e19k9xsayI/BwRMyNiLeAK4D9y8rsD1yU318N7JX/Yt4fuCIi3oyIecDDub5a6uyUNkfEtIh4Km+fDfSRtFod21b3NlerUNKGwFoRcW+k/0MvBoY0YJuH52M7QqttjohHI2ImsLjs2E8Ct0bEixHxEnAr8KlGuM7V2hwR/46Ih/L7p4BngYqZVepsRa5zRfn3Zk/S7xGk36uGuM5lhgJ/i4jX69g2B8YG0x94vPD5ibytYpmIeAd4BXhvC8fWUmdntbnoQOC+iHizsO3P+XbIj+p8u2xF2/wBSdMk3SFpt0L5J1qpszPbXDIMuLxsW2de57Ye2wjXuVWSdiT1hB4pbP5pvsV6Zp3/AFzRNveWNEXSvaVbkqTfm5fz79Hy1Nmaev27dDDL/j6v8HV2YLROJ2kb4OfA/xY2j4iIjwC75deXOqNtFTwNbBIR2wHfBi6TtFYnt6kmkj4GvB4R9xc2N+p17rJyr/YS4MsRUertfA/YEtiBdPvvhE5qXiWbRkq19kXgN5I27+wG1SJf548A4wqb63KdHRgby5PA+wufN87bKpaRtDKwNvBCC8fWUmdntRlJGwPXAYdGxLt/XUfEk/nna8BlpFsvnd7mfKv6hdy2qaQewYdy+Y1bqbNT2lzYv8xf1w1wndt6bCNc56ryH0l/BX4QEfeWtkfE05G8CfyZxrnOxd+BuaRnztuRfm/65d+jNtdZg3r8u/QF4LqIeLu0oV7X2YGxsUwGtsijwVYl/UM2tqzMWKA0Qm8oMD4/axkLHKw0MvEDwBakQQq11NkpbZbUj/SPyIkRcXepsKSVJa2b368CfAa4n/pZkTavJ6lXbttmpOs8NyKeBl6V9PF8O/JQ4C+N0Obc1pVI/5C8+3yxQa5zNeOAfSW9R9J7gH2BcQ1ynSvK5a8DLo6Iq8v2bZh/ivSsriGuc76+q+X36wK7AP/Kvze3k36PIP1eNcR1LhhO2R96dbvOKzp6x6/6voD9gH+TeiI/yNtOAT6X3/cGriINrpkEbFY49gf5uDkURupVqrMR2gz8EFgATC+81gfWAKYCM0mDcn4L9GqQNh+Y2zQduA/4bKHOQfl/xEeAc8iZpTq7zXnfYODesvoa4TrvQHq+tIDUS5ldOPbw/F0eJt2WbJTrXLHNwCHA22W/zwPzvvHArNzuS4G+DdLmnXO7ZuSfXynUuVn+PXo4/16t1ghtzvuaSD3MlcrqrMt1dko4MzOzAt9KNTMzK3BgNDMzK3BgNDMzK3BgNDMzK3BgNDMzK3BgNGtQWrLqwf2SbsjzPls7Zn4r+/tJ+lrh80aSrm7pmBrb2iSpnnPzajnnQEn7deQ5rWdwYDRrXAsjYmBEDCCtIHB0HersR1qFA0jJriNiaAvlG1LOyDKQNBfOrK4cGM26hokUkixLOl7S5Jws+eTywpL6SrpNaT27WZJKKxecDmyee6JnFHt6OYn0NoU6JkgaJGkNSRdImpSTp7e4OoukkZKul3SrpEclfV3St/Ox90pap1D/bwu94h3z9nXy8TNz+W3z9pMkXSLpblIu0lOAYfn4YZJ2VFoPcZqkeyR9uNCeayXdrLS24y8Kbf1UvkYzJN2Wt7Xp+1o3VM9MBn755Vf9XuR18khr110FfCp/3hcYRVqPcCXgRmD3smNWJi3PBLAuKXuJKFtPsvgZ+BZwcn6/ITAnv/8ZcEh+34+UrWSNsrYW6xmZz7cmadmlV4Cj8r4zgWPz+wnA+fn97oXjzwZ+nN/vCUzP708iZerpUzjPOYU2rAWsnN/vDVxTKDeXlDu2N/AYKU/neqQVHj6Qy61T6/f1q3u/Sglizazx9JE0ndRTfIC0JiGkwLgvMC1/7kvK2Xpn4VgBP5O0O2k9u/7ABq2c70rgFuDHpLyqpWeP+wKfk3Rc/twb2CS3qZrbIyUmf03SK8ANefss0uKzJZcDRMSdktbKz1F3JaXeIyLGS3qvlqxgMjYiFlY559rARZK2AAJYpbDvtoh4BUDSv4BNSQsg3xlp/VIi4sUV+L7WjTgwmjWuhRExUNLqpKTaRwNnkYLeaRFxXgvHjiD1iJoj4m1Jj5L+ga8qIp6U9EK+dTkMOCrvEnBgRMxpQ9uL62ouLnxezNL/7pTnpGwtR+WCFvadSgrIB0hqIvVIK7VnES3/27c839e6ET9jNGtwkVYnPwb4Th50Mg44XFJfAEn9Ja1fdtjawLM5KO5B6iEBvEa6xVnNGOC7wNqRVk8nn+8becUCJG1Xj++VDct17gq8knt1d5ECO5IGA89HxKsVji3/LmuzZOmikTWc+15gd6XVaCg9+6R9v691AQ6MZl1AREwjrYIxPCJuIa2dOFHSLNItz/JgNxoYlPcfCjyY63kBuDsPdjmjwqmuJi0BdGVh26mk25IzJc3On+vlDUnTgHOBr+RtJwHNkmaSBgsdVuXY24GtS4NvgF8Ap+X6Wr0bFhHPAUcC10qaQfqjANr3+1oX4NU1zKxTSJoAHBcRUzq7LWZF7jGamZkVuMdoZmZW4B6jmZlZgQOjmZlZgQOjmZlZgQOjmZlZgQOjmZlZwf8H9OeqZXvGn5oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.9702253302253302"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf = RandomForestClassifier(max_depth=16, random_state=0, n_estimators=10)\n",
    "rf.fit(X_train, Y_train)  \n",
    "print(rf.feature_importances_)\n",
    "importances = rf.feature_importances_\n",
    "indices = np.argsort(importances)\n",
    "features = X_train.columns\n",
    "plt.title('Feature Importances')\n",
    "plt.barh(range(len(indices)), importances[indices], color='b', align='center')\n",
    "plt.yticks(range(len(indices)), [features[i] for i in indices])\n",
    "plt.xlabel('Relative Importance')\n",
    "plt.show()\n",
    "sklearn.metrics.accuracy_score(Y_test, rf.predict(X_test))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.array(X_train.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Greenwell et al.: \"A Simple and Effective Model-Based Variable Importance Measure\"\n",
    "import math\n",
    "from sklearn.inspection import partial_dependence\n",
    "\n",
    "# variable feature importance\n",
    "def vfi(model, data, a, k):\n",
    "    pd, _ = partial_dependence(model, data, [a])\n",
    "    imp = 0\n",
    "    if pd.shape[1] < 100 and pd.shape[1] < data.shape[1]: # categorical\n",
    "        imp = (max(pd[0]) - min(pd[0]))/4\n",
    "    else: # continuous\n",
    "        for i in range(k):\n",
    "            right = 0\n",
    "            for j in range(k):\n",
    "                right += pd[0][j]\n",
    "            imp += (pd[0][i] - right/k)**2\n",
    "        imp = math.sqrt(imp/(k-1))\n",
    "    return imp\n",
    "\n",
    "# standard deviation of feature importance values of feature b given a certain value of feature a\n",
    "def std(model, data, b, a):\n",
    "    unique = np.unique(data[:,a])\n",
    "    iv_ba = []\n",
    "    for uv in unique:\n",
    "        uv_a = data[np.where(data[:,a] == uv)]\n",
    "        iv_ba.append(vfi(model, uv_a, b, 4))\n",
    "    iv_ba = np.array(iv_ba)\n",
    "    return np.std(iv_ba)\n",
    "    \n",
    "# feature interaction based on the median of the standard deviation of conditional feature importances between two features\n",
    "def fint(model, data, a, b):\n",
    "    sb = std(model, data, b, a)\n",
    "    sa = std(model, data, a, b)\n",
    "    return (sa+sb)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc4AAAEWCAYAAADvi3fyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debxd0/3/8ddbzGKooWooqVBD0ogkomZRtHx9K4iiKdL2S7VK9VeUjqoDbakWNbYEDUIMRdVQxKyRyIyYQqkpxiZpYkg+vz/WOrJznHPvOXc8uff9fDzOI+fsvfba6+zr+ty191qfpYjAzMzMarNUZzfAzMxsSeLAaWZmVgcHTjMzszo4cJqZmdXBgdPMzKwODpxmZmZ1cOA0ayOSNpA0R1KPGsruIunFJvaPlPSLtm2hmbUFB07rliTdKumUCtv3kfSKpKXrrTMi/hURPSNiQdu0smUkhaS5OYjPkfR2G9Q5QtL9bdG+Os45VtL/deQ5q/EfMlbkwGnd1aXAVySpbPshwKiI+KCeyloSaNvZljmI94yI1Tq7MQ14fWpWyx0E614cOK27ugFYA9ixtEHSx4C9gcvy5/+RNFHSfyS9IOnkQtleuWf3dUn/Au4qbFs6l/mqpMclzZb0rKRvlDdC0g8kvS7pOUnDqzVW0t6SJkl6W9KDkvq15Es3VY+kEyU9k9v7mKR98/bNgfOBbYs92PIeYXmvNF+LoyQ9BTyVt20m6Q5Jb0qaIelLNbZ7F0kvSjpB0muSXpY0VNJekp7M9f2gUP5kSWMkjc7f51FJWxb2b57b/7ak6ZK+WNg3UtJ5km6RNBf4OjAcOCF//5uaul7FayHpdElvSZopac/C/tUlXSLppbz/hhp/Rt+X9O98zhmSPlfL9bM2FhF++dUtX8BFwJ8Kn78BTCp83gX4DOkPzH7Aq8DQvK8XEKQguxKwQmHb0rnM/wC9AQE7A/8FBhTq/gD4HbBc3j8X2DTvHwn8Ir/fCngN2AboARwGPAcsV+V7BbBxhe1N1gMcAKybv++BuT3r5H0jgPvL6hsL/F/h82JlcjvuAFbP12cl4AXgq8DSuT2vA1tU+R4f1l+4Xj8BlgEOB2YBVwArA32AecCncvmTgfeBYbn8ccDM/H4Z4GngB8CywK7A7LJr/w6wfb4Wyxd/HoX2NXe93s/t7AF8E3gJUN7/N2A08LHcnp2b+xkBm+brt27hv8Henf171B1f7nFad3YpMEzS8vnzoXkbABExNiKmRsTCiJgCXEkKcEUnR8TciJhXXnlE/C0inonkHuB2Cj3c7McR8W7e/zegUg/sCOCCiPhnRCyIiEuBd4HPNvHdHs09lrclnVVLPRFxTUS8lL/vaFIvcXAT56jFqRHxZr4+ewPPRcQlEfFBREwEriUFoFq8D/wyIt4HrgLWBP4QEbMjYjrwGLBlofyEiBiTy/+OFAA/m189gdMi4r2IuAu4GTi4cOxfI+KBfC3mV2pMDdfr+Yi4KNIz70uBdYC1Ja0D7AkcGRFvRcT7+ecPTf+MFpAC6BaSlomI5yLimRqvnbUhB07rtiLiflKPZ6ik3qT/6V1R2i9pG0l3S5ol6R3gSNL/rIteqFa/pD0lPZxvI74N7FV2/FsRMbfw+XlSD6bchsD3CoHwbeCTVcqWDIiI1fLrmFrqkXRo4Rbh20DfCt+3XsXrsyGwTdn5hwOfqLGuN2LRwKvSHyqvFvbPIwXEj5w7IhYCL5K+67rAC3lbyfPAelXaXVEN1+uVwvn/m9/2JF3zNyPirQrVVv0ZRcTTwLGk3vRrkq6S1NR/A9ZOHDitu7uM1NP8CnBbRBT/R3wFcCPwyYhYlfScr3wwUcXlhSQtR+pNnQ6sHWmAzi1lx39M0kqFzxuQbueVe4HU01qt8FoxIq6s+Vs2U4+kDUm3rr8NrJHbO63Q3krfcy6wYuFzpQBYPO4F4J6y8/eMiG/W+T1q9cnSG0lLAeuTru9LwCfztpINgH9XafdHPtdwvZryArC6pEqDtpr8WUfEFRGxAynABvDrGs5nbcyB07q7y4DdSM+iLi3btzKpZzBf0mDgy3XUuyzpttos4IM8MGSPCuV+JmlZSTuSbmVeU6HMRcCRuQcsSSspDVxauY72NFfPSqT/Ec+CNLCJ1IMqeRVYX9KyhW2TgP0krShpY9IgmqbcDHxa0iGSlsmvrZUGH7WHgZL2UxqsdSzplufDwD9Jz5tPyG3YBfhf0u3fal4FNip8bu56VRURLwN/B86V9LHchp3y7qo/I0mbSto1/1E2n9TDXljlNNaOHDitW4uI54AHSf8jvLFs97eAUyTNJg1KubqOemcDx+Rj3iIF3fL6X8n7XgJGkZ55PVGhrvGkwH5OLv80afBJXZqqJyIeA84AHiIFic8ADxQOvwuYDrwi6fW87UzgvVz+0vwdmjr/bNIfDweRvvMrpB7TcvV+lxr9lTRo5y3SNKP98vPE90iBck/SrfpzgUMrXfuCP5OeLb4t6YYarldzDiE9s32CNBjoWGj2Z70ccFpu8yvAx4GT6jintZHSCC8zsy5DaerQxhHxlc5ui3U97nGamZnVwYHTzMysDr5Va2ZmVgf3OM3MzOqwxCZetkXWXHPN6NWrV2c3w8xsiTFhwoTXI2KtlhzrwNkF9OrVi/Hjx3d2M8zMlhiSnm/psb5Va2ZmVgcHTjMzszo4cJqZmdXBgdPMzKwODpxmZmZ1cOA0MzOrgwOnmZlZHRw4zczM6uAECF3AhAmgWtadNzPrIjozzbp7nGZmZnXoFoFT0g8lTZc0RdIkSdtIGitpRv78uKQjyo7pLykkfaHGcwzN5Tdrn29hZmaNoMsHTknbAnsDAyKiH7Ab8ELePTwi+gPbA7+WtGzh0IOB+/O/tai3vJmZLYG6fOAE1gFej4h3ASLi9Yh4qaxMT2AusABAkoADgBHA7pKWb+oEknoCOwBfBw4qbN9F0r2S/pZ7t+dLWirvmyPpzNwTvlPSWpI2kzSucHwvSVOrnPMISeMljYdZdV0QMzNrue4QOG8HPinpSUnnStq5sG+UpCnADODnEbEgb98OmBkRzwBjgf9p5hz7ALdGxJPAG5IGFvYNBo4GtgB6A/vl7SsB4yOiD3AP8NOIeAJYVtKncpkDgdGVThgRF0bEoIgYBC1aGcfMzFqgywfOiJgDDASOIHXNRksakXcPz7dvNwCOk7Rh3n4wcFV+fxXN335tqvy4iHg2B+UrST1TgIUsCop/KWy/mhQwoYnAaWZmnaNbTEfJQWssMDbf+jysbP8sSY8C20h6Edgf2EfSDwEBa0haOSJml9ctaXVgV+AzkgLoAYSk40vVlzenWjPzv6OBayRdl5oWT9X5dc3MrB11+cApaVNgYSEA9QeeB/oWyqwIbAX8BvgcMCUiPl/YfymwL3BZhVMMAy6PiG8Uyt8D7Jg/Ds63Xp8n9SAvzNuXysdeBXyZNLCIiHhG0gLgx9TY2xw4ELyOtZlZx+jyt2pJA38ulfRYfp65BXBy3jdK0iRgAjAyIiaQbrNeX1bHtVS/Xdtc+UeAc4DHgZmFsnNJQXUaqcd6SuH40cBXSLdtzcysgSg6M/1CFydpF+C4iNi7wr45EdGzbc4zKMBdTjNrfI0SciRNSIMr67dE9DglfULSVZKekTRB0i15OsbNLahraE6E8LikqZKGFvZtlhMiTJTUW9IxudwoSSMkzcr7H5N0eJX6B0k6q7am6ApJz+bv9JCkfev9PmZm1rEaPnDmOZXXA2MjondEDAROAtZuQV1bAqcD+0TE5sAXgdMl9ctFhgJjImKrPBXlW8DuETGcNH2kZCFwnqRpktYo1L90RIyPiGMAImJsld6mgCnAvRGxUf5OBwHrVyjb5Z9Dm5ktSZaE/ykPAd6PiPNLGyJisqSPAZ+TNIY00GcC8JWICEl7Ab8jPUd8ANgoB7DjgF9FxMxcz0xJpwLHS7oSOBZYIOlzpLmdGwF/l3Qx8BYwOiK+DSDpYeAY4GhJvXPZf0m6gHx7VtLJpKkuG+V/fx8RZ5Geab5X9p2eB87OdY8gzffsSRqlW5x7amZmnWhJCJyloFjJVkAf4CVSgNw+ZdLhAmCnHBivLJTvQ+pxFo0HjoqIWySdD8yJiNMBcp7aIRHxemHuJ5I2IgXDp/OmLYAdImJefq5ZtBkp+K8MzJB0Xm7Ho8187wFAv4h4s9LOnFs359fdoJmqzMysrTT8rdpmjIuIFyNiITAJ6EUKVM+WepWkpANt5cA8CvdK4BuFoHZjRMyrcszfIuLdiHgdeI0Kt5gl/VHSZEmPFDbfUS1ogjMHmZl1liUhcE4nZf6p5N3C+wU034N+rEJdA/M5ajE6IvpHxDYRUZyCMreJYyq1cTqpRwlARBxFmj9ajIBN1WlmZp1kSQicdwHLqbDsVx7Ms2OV8jOAjST1yp8PLOw7HTiptC//+wPgjDZsby3uApaX9M3CthU7uA1mZtYCDf+MMw/22Rf4vaTvA/OB54AbqpSfJ+lbwK2S5pISEJT2Tcp13CRpGeB94ISImNTe36OsjZGnwZwp6QRSDt25wPdbUp8zB5mZdZwumQBBUs+ImJOnffwReCoizuzsdrWXQYMGxXhHTjOzmrUmAULD9zhb6HBJhwHLAhNJo2y7rAkTQOrsViy5uuDfjmbWjpaEZ5x1i4gz8yCeLSJieET8t946JK1dT2YfScdKmi9p1da13szMGlmXDJytlW/x3kAzmX3KsvocTHqeuh9mZtZlOXBWVjGzT0ScnXPW3ijpLuBOgJw5qCfwIwqrqOSyf5U0VtJTkn6at/eS9ETOgfu4pDGSVpT0BUnXFI7fpSX5eM3MrP04cFbWXGafAcCwiCilwjuItK7mfcCmkopJDgaTFsbuBxwgqfQwelPg3Jwz9z+kvLj/IC2mXcqLe2Cu9yNykvvxKVPSrLq/oJmZtYwDZw0qZPYpz+pzMHBVzmB0LXBAYd8dEfFGzix0HbBD3v5CRDyQ3/+FlLLvA+BW4H/zbeD/Af5aqU3OHGRm1jm66qja1ppO6iUCKbOPpDVZtOjlh1l9JH0G2AS4Iz0aZVnSgtXnlA4vqzua2X4V8G3gTWB8RMxu1TcxM7M25R5nZfVk9jkYODkieuXXusC6kjbM+3eXtLqkFUjLlpV6mRtI2ja//zJwf35/D+lW8OFUuU1rZmadx4GzgkhZIYYCO0uaKWkccCmVM/scRFovtOj6vB1gHOn27RTg2ogo9VpnAEdJehz4GHBePvcC4GZgz/xvswYOTHMR/WrZy8ysHl0yc1CjyEuRDSqt4VnY3gu4OSL6tsV5nDnIzKw+zhzUzXWVzEH+G87MlgTN3qqVtEDSJEnT88jS70laKu8bJOmsJo7tMvMQJf2g7HOzmYUiYmR5bzNvf66teptmZtaxannGOS+nr+sD7E569vZTgIgYHxHHtFfjyjLzdLYPA2etmYVy2Ub6DmZm1kp1DQ6KiNeAI4BvK/mwRylp59wznSRpoqSV82GrSPqbpBmSzi/0VueU6pU0TNLI/H5kLvdP4DfV6pV0vKRHJE2R9LO8rZSRZ6SkJ3Nmnt0kPZAz9wzO5VaSdLGkcbnOffL2EZKuk3RrLv+bvP00YIXchlE0kVmoUM+H2YXytfqtpGmSpko6MJdbR9K9ud5pknaU1CO3v1T2u/X8jMzMrH3V3RuKiGcl9QA+XrbrOOCoiHhAUk/SupmQMudsATxPmty/HzCmmdOsD2wXEQsk3VRer6Q9SHMnBwMCbpS0E/AvYGNSAoKvkXLHfpmUdOCLpF7jUOCHwF0R8TVJqwHjJP0jn7s/sBXwLjBD0tkRcaKkb0dEfwBJx9B0ZiFIU0r6RcSbkvbP9W4JrAk8Iune3LbbIuKX+ZqumMutV7qVm9v3EUoLe+fFvTdopilmZtZW2nI6ygPA73JQWS1nwQEYFxHP5mkWV7Ioc05Trsnlq9W7R35NJAWwzUiBFGBmREzNWXymA3fm6SVTgV65zB7AiZImAWOB5VkUfe6MiHciYj7wGFCaj1mVPppZCBbPLrQDcGVELIiIV0lzNbcmBfavSjoZ+ExOdvAssJGksyV9gZSO7yOcOcjMrHPUHTglbQQsAF4rbo+I04D/A1YAHpC0WWlXWRWVMucsX1bmw8w8VeoVcGp+9to/IjaOiD/nQ94t1LOw8Hkhi3rYAvYvHL9BRDxe4fgFVO6VTyf1KEttPAr4HItHsLnlB5WLiHuBnYB/AyMlHRoRb5F6pmOBI4E/NVePmZl1nLoCp6S1gPOBc6JsAqik3rmn92tST6oUOAdL+lR+tnkgizLkvCpp87y9qXUuK9V7G/C1fOsWSetJKr913JTbgKPzIB8kbVXDMe9LWia/ryezEKTk7wfm55drkYLlOKXsQq9GxEWkADlAKbXfUhFxLWm1lQFVazUzsw5XyzPOFfItzWWAD4DLgd9VKHespCGknt104O/AtqRgdw7p2ePdLMqycyIpM84sUg7YnlXO/5F6I+JdSZsDD+XYNwf4CqmHWIufA78HpuTAPRPYu5ljLszlH42I4ZKGAmdKOiF/h7lUziwE6TtvC0wm9bRPiIhXJB0GHC/p/fwdDgXWAy4pDaICTmruywwcCM5/YGbWMZw5qAtw5iAzs/rImYO6t0bIHOS/v8ysu3CS9zYg6YdKmZWm5DmZ20g6VlJTzz1Lx9ZUzszMGoMDZyspLQ22NzAgIvoBuwEvAMfS9IChklrLmZlZA3DgbL11gNcj4l2AiHgdGAasC9wt6W4ASedJGp97pqVMR8dUKLeHUt7bRyVdUxo5bGZmjcGDg1opB7b7Sb3GfwCjI+IeSc+RlhR7PZdbPWcR6gHcCRwTEVOK5fJUlOuAPSNirqTvA8tFxCkVzlvMHDQwJWbqPP7PyMyWJB4c1IkiYo6kgcCOwBBgtKQTKxT9Ug52S5N6qVuQFrcu+mze/kCeZrMs8FCV815ImiKDNMhhy8ysgzhwtoGcHnAsMFbSVOCw4n5JnyLl8t06It5SSmhfni0JUkajOyLi4PZtsZmZtZSfcbaSpE0lbVLY1J9033Q28OEKMaQECe9IWpu0NFtJsdzDwPaSNs51ryTp0+3ZfjMzq497nK3XEzg7r2LyAfA06dnjwcCtkl6KiCGSJgJPkEbcPlA4/sKyciOAKyUtl/f/CHiyqQY4c5CZWcfx4KAuwJmDzMzq48FB3VxnZg7y311m1t10y2eckkLSGYXPx+U1MVtT5y6S3smZg0qv3Zoov66kqgt6S1pN0rda0yYzM2t73TJwktbc3C/Pm2xL9xXW+OwfEf+oVjAiXoqIYU3UtRrgwGlm1mC6a+D8gDQo57vlOyStJelaSY/k1/Z5+9TcC5SkNyQdmrdfJmn3aieStHXOYbt8HiU7XVJfSb0kTctl+kgal3upU/Io3dOA3nnbb9vjIpiZWf268zPOP5LW1/xN2fY/AGdGxP2SNiAter05aSTs9qSpJs+SEh5cRlpn85vA1sCOee3Skv0j4hFJNwK/AFYA/hIR0yT1KpQ7EvhDRIyStCzQg7Read+I6F+p8WWZg1p0AczMrH7dNnBGxH8kXQYcA8wr7NoN2EKLRtusktPq3QfsRAqc5wFHSFoPeCunx4N0q7bSgtinkBb0np/PV+4h4IeS1geui4in1MxoH2cOMjPrHN31Vm3J74GvAysVti0FfLbwnHK9iJgD3EvqZe5IyhI0i5TM/b4azrMGab7nylTIGBQRVwBfJAXwWyTt2uJvZGZm7apbB86IeBO4mhQ8S24Hji59kNQ/l30BWBPYJCKeJSV2P44UUJtzAfBjYBTw6/KdkjYCno2Is4C/Av1YPKOQmZk1iG4dOLMzSAGx5BhgUB6k8xjp+WPJP1mUxec+YD1SAC3ZsWw6yrA8iOj93Ks8Ddi6Qo/yS8C0/Hy0L3BZRLxBSvY+rbnBQQMHpvmUnfEyM+tunDmoC3DmIDOz+jhzUDfXnpmD/HeVmdniuvStWkkL8i3TyZIelbRdK+oaK+kjf51I+lqe4zkl31bdJ28fIWndGuqtqZyZmTWGrt7jnFeaBynp88CpwM5tVXmePvJDYEBEvJOnrayVd48ApgEvNVNNreXMzKwBdOkeZ5lVgLcAJPWUdGfuhU4t9BJ7SXpc0kU5w8/tklYoViJpKUkjJf0C+Dhp9OscgIiYExEzJQ0DBgGjco93BUk/yZmIpkm6MGcgqlRuoKR7JE2QdJukdTruEpmZWXO6euBcIQekJ4A/AT/P2+cD+0bEAGAIcIYWZRzYBPhjRPQB3gb2L9S3NGlKyVMR8SNgMvAqMFPSJZL+FyAixgDjgeF5Lug84JyI2Doi+pIyCO1dXo6UCvBsYFhEDAQuBn5Z6YtJOkLSeEnj05RSMzPrCN3pVu22wGWS+gICfiVpJ2AhaVrJ2vmYmRFRSps3AehVqO8C4OqI+CVARCyQ9AVSur3PAWdKGhgRJ1doyxBJJwArAqsD04GbyspsSpqOckeO4z2Alyt9MWcOMjPrHF09cH4oIh5SWg1lLWCv/O/AiHhf0nMsyujzbuGwBaTeYcmDpAB4RkTMz/UGMA4YJ+kO4BLg5OK5JS0PnAsMiogXlJYw+0gGIVJAnx4R27bmu5qZWfvp6rdqPyRpM1IP7g1gVeC1HDSHABvWWM2fgVuAqyUtrbSm5oDC/v6kXLaweOafUpB8PQ8gKi4nViw3A1gr946RtIykPjV/STMza3ddvce5QmG1EgGH5duro4CbJE0lPWN8otYKI+J3klYFLietYHJ6nk4yn/SwsZRpaCRwvqR5pBVULiKNnn2FlPCdKuWGAWflcyxNyqc7vak2DRwIzn9gZtYxnDmoC3DmIDOz+jhzUDfXksxB/nvJzKxlus0zzs4mad+yBPCTJC2UNFzSmM5un5mZ1ca3ajuJpCOA4cCQiFjYuroGRXpUWzv/2M2sO2vNrVr3ODuBpE8DPwEOATaQNC1vHyHprzkv7lOSftqpDTUzs49w4OxgkpYBrgC+FxH/qlBkMClbUT/ggEqJ5XM9zhxkZtYJHDg73s9JSQ5GV9l/R0S8kdP0XQfsUKlQRFwYEYPSrYa1KhUxM7N24FG1HUjSLqTe5IAmipU/ffTTSDOzBuIeZweR9DFSOr5DI2J2E0V3l7R6XpVlKPBAhzTQzMxq4h5nxzmStAzZeVp80uWVZeXGAdcC6wN/iYhmh8s6c5CZWcdx4OwgEXEqaSHtSn5deP9iRAztgCaZmVkLOHB2AfVmDvIcTjOzlusSzzglfULSVZKekTRB0i15usbNLahrqKQpkh6XNFXS0MK+zXLGn4mSeks6Jpcbledgzsr7npJ0m6TtajzfFgARMTIivl1vm83MrOMs8YFT6YHh9cDYiOgdEQOBk1i0MHU9dW0JnA7sExGbA18krX7SLxcZCoyJiK0i4hngW8DuETE87x+d920CnAZcJ2nzZk47FNii3raamVnnWOIDJzAEeD8izi9tiIjJwH1AT0ljJD2Re4UCkLRX3jZB0lmFnulxwK8iYmauZybpueTxkvYCjgW+KeluSecDGwF/l/Td8kZFxN3AhcAR+ZyHS3pE0mRJ10paMfdIvwj8Nvdke+fXrblt9+V1RM3MrEF0hcDZF5hQZd9WpGC3BSnIbS9peeACYM/cOy1mD+hToa7xQJ+IuAU4HzgzIoZExJHAS6Rcs2dWOf+jQCnwXRcRW0fElsDjwNcj4kHgRuD4iOife7EXAkfnth0HnFupYmcOMjPrHF19cNC4iHgRIC9o3QuYAzxb6lWSpoMc0U7nLw7Z6SvpF8BqQE/gto8UlnoC2wHXFKasLFep4oi4kBRkc5J3MzPrCF0hcE4HhlXZ927h/QKa/76PAQOByYVtA/M5WmIrUu8SYCQwNCImSxoB7FKh/FLA2xHRv4XnMzOzdtYVbtXeBSyXl+kCIA/m2bFK+RnARpJ65c8HFvadDpxU2pf//QFwRr2NkrQzqSd7Ud60MvByTvI+vFB0dt5HRPwHmCnpgFyH8oAlMzNrEEt84Iy0oOi+wG55Osp00oCeV6qUn0caDXurpAmkwPVO3jcJ+D5wk6QngJuAE/L2WhyYB/k8SQq4+0dEqcf5Y+CfpBR6TxSOuYo0+GiipN6koPp1SZNJPd19mjvpwIFpbmatLzMza7luuZC1pJ4RMSePsv0j8FQTA3wa3qBBg2K8c+6ZmdXMC1nX7/A8WGg6sCpplO0Sq5Q5qNLLzMzaVrcMnBFxZp7+sUVEDI+I/7amPkkL8i3ayZIerSVjUBN1ja22eLWZmXW+rjCqthHMK42ElfR50jPWnTu3SWZm1h66ZY+zna0CvAXpWaqkO3MvdKqkffL2XjnH7UWSpku6Pa+/+SFJS0kamed+mplZg3CPs22skJ+ZLg+sA+yat88H9o2I/0haE3hY0o153ybAwRFxuKSrgf2Bv+R9SwOjgGkR8ctKJ8zTb/IUnA3a/huZmVlF7nG2jXn5melmwBeAy/KIXQG/kjQF+AewHouSz88sTHOZQMpqVHIBTQRNSJmDImJQGhW2VrViZmbWxhw421hEPASsSYpmw/O/A/Mz0FdJvVJoOqvRg8CQnFfXzMwaiANnG8urmfQA3iBNdXktIt6XNATYsMZq/gzcAlwtybfTzcwaiANn21ghT0eZBIwGDouIBaTnlIMkTQUOZfGMQU2KiN8BE4HLJTX5c2oqc5CZmbWtbpk5qKtx5iAzs/q0JnOQbwN2AaXMQZX47yIzs7blW7VNkBSSzih8Pk7Syc0cc7Kkf5du3ebXau3eWDMz6xAOnE17F9gvz8GsRymlX+n1dns0zszMOp4DZ9M+AC4Evlu+I2f/uUvSlJwdqMksBJJGSLpB0h2SnpP0bUn/Ly8n9rCk1XO5sZL+kHuq0yQNbp+vZmZmLeHA2bw/AsMlrVq2/Wzg0ojoRxo9e1Zh33cLt2nvLmzvC+wHbA38EvhvRGwFPEQadVuyYp73+S3g4kqNknSEpPGSxsOs1nw/MzOrgwNnMyLiP8BlwDFlu7YFrsjvLwd2KOwr3qodUth+d0TMjohZpMWzb8rbp7J45qAr87nvBVap9IzUmYPMzDqHA2dtfg98HViplfUUswUtLHxeyHr5IQ0AABsBSURBVOIjnMvHwnpsrJlZg3DgrEFEvAlcTQqeJQ8CB+X3w4H72vCUBwJI2gF4JyLeacO6zcysFRw4a3cGKQdtydHAV3MC90OA7xT2FZ9xTpLUq85zzZc0ETifxYN1Rc4cZGbWcZw5qMFIGgscFxE1pwJy5iAzs/o4c1A3VylzkP8eMjNrH75V284kfULSVZKekTRB0i2SdpI0Ju/vL2mvUvmI2KWe3qaZmXUsB852lBezvh4YGxG9I2IgcBIQETEsF+sP7FWtDjMzaywOnO1rCPB+RJxf2hARk4EXclagZYFTgAPzIKIDJT0laS0ASUtJerr02czMOp8DZ/vqC0yotjMi3gN+AozOyRJGA38hTW8B2A2YnBMmLMaZg8zMOocDZ+O5mEXp974GXFKpkDMHmZl1DgfO9jUdGFjPARHxAvCqpF2BwcDf26NhZmbWMg6c7esuYDlJR5Q2SOoHfLJQZjawctlxfyLdsr0mIha0eyvNzKxmDpztKFJ2iX2B3fJ0lOnAqcArhWJ3A1uUBgflbTcCPalym7ZcpcxBZmbWPpwAoZ1FxEvAlyrs6pv3v0laZqxoS9KgoCfauXlmZlYnB84GI+lE4JssGlnbLGcOMjPrOL5VW0bSnDasa4SkWZIm5vmZt0narqljIuK0iNgwIu5vq3aYmVnbceBsf6MjYquI2AQ4DbhO0uad3SgzM2sZB84aSOol6S5JUyTdKWkDST0kzVSymqQFknbK5e+VtEl5PRFxN3AhcEQu11/Sw7ne6yV9LG/fWNI/JE2W9Kik3h35fc3MrDoHztqcDVwaEf2AUcBZeZrIDGALYAfgUWBHScsBn4yIp6rU9SiwWX5/GfD9XO9U4Kd5+yjgjxGxJbAd8HJ5Jc4cZGbWORw4a7MtcEV+fzkpUALcB+yUX6fm7VsDjzRRlwAkrQqsFhH35O2XAjtJWhlYLyKuB4iI+RHx3/JKnDnIzKxzOHC2zr3AjqQMP7cAqwG7kAJqNVsBj7d7y8zMrF04cNbmQeCg/H44iwLjONKt1IURMR+YBHyDFFA/QtLOpOebF0XEO8BbknbMuw8B7omI2cCLkobmY5aTtGI7fCczM2sBB86PWlHSi4XX/wOOBr4qaQopwH0HICLeBV4AHs7H3kdKnze1UF9pybAngR8A+0dEqcd5GPDbXG9/0hJj5HMck7c/CHyiqQY7c5CZWcdR+P+yS7xBgwbF+PHjO7sZZmZLDEkT0hiR+rnH2QWUMgeVZw8yM7O21+UDZ55fOUnSNEk3SVqthfX0kjStwvZTJO3WivYNzvM+Z+QMQ3/yM00zs8bV5QMnMC8i+kdEX+BN4Ki2rDwifhIR/2jJsZLWBq4hzeXcNCK2Am6lbJkxSc4pbGbWILpD4Cx6CFiv9EHS8ZIeyZl7fpa3nSbpqEKZkyUdV61CSSMlDcvvn5P0s5ztZ6qkzfL2lSRdLGlc7lXukw8/ipRY4aFSfRExJiJezee9XNIDpLmjZmbWALpN4JTUA/gcaa1LJO0BbEKag9kfGJhT5o1m8WXAvpS31er1iBgAnAeUAu4PgbsiYjAwhDSSdiXS0mITmqhrC2C3iDi4wvdx5iAzs07QHQLnCpImkRaPXhu4I2/fI78msigN3iYRMRH4uKR1JW0JvBURL9RxvuvyvxOAXoVznZjbMRZYHtighrpujIh5lXY4c5CZWefoDs/O5kVE/zzg5jbS7dGzSKnvTo2ICyoccw0wjDR/sp7eJsC7+d8FLLq+Is3fnFEsKGk6MBD4a5W65tZ5bjMza2fdoccJQM73egzwvTzY5jbga5J6AkhaT9LHc/HRpExBw0hBtLVuA46WVMpTu1Xefg5wmKRtSgUl7ZcHDZmZWQPqNoETIN+GnQIcHBG3kxK3PyRpKjCGPJo1Iqbn9/+OiOLKJJuWZRU6oMZT/xxYBpiSe5k/z+d5lRSgT8/TUR4HPg/Mrud7FTMHmZlZ+3LmoC7AmYPMzOrjzEHdXClzkJmZtb8uEzgrZfapNgezOPeyDc47Nt9mnZznhPZvi3pz3btIurmt6jMzs9brMoGzVu2UhWd4RGwJnAv8tgHaY2Zm7aRbBM7cK/x9ShaQlgQDdssJBJ6UtHcu10vSfTnzz6OStsvbd8l1jJH0hKRRpRGyZcozE80pvB8maWR+P1LS+ZL+Cfwm56t9KGcVelDSpu1zJczMrLW6U29n2dKD4BzAepGyBvUG7pa0MfAasHtEzJe0CXAlUHp4vBXQB3gJeADYHri/7BxfAG6osT3rA9tFxAJJqwA7RsQHOWH8r4D9mzpY0hGkRbGpLZeCmZm1ha4UOKsNDy5tL09kcHVELASekvQsKXPQTOCc/JxyAfDpQvlxEfEiQM4A1ItFgXOUpGWBnqT0fbW4JiIW5PerApfmYB2kqStNiogLgQtTewZ5aLSZWQfpSrdq3wA+VrZtdeD1/L48C095sAngu8CrwJaknuayhf3vFt4XswIBDAc2Ai4Fzq5yjuXLzldsz8+Bu/MKLv9boayZmTWILhM4I2IO8LKkXQEkrU66dVp+O7XkAElLSepNCnozSD2/l3NP9BCgRx3nD+DHwGdLq6IAr0raXNJSwL5NHL4q8O/8fkSt5zQzs47XZQJndijw43wr9S7gZxHxTJWy/wLGAX8HjoyI+aRRsYdJmky6dVtXrtickP0M4Pi86UTgZuBB4OVqxwG/AU6VNJEW3D4vZQ4yM7P258xBXYAzB5mZ1ceZg7o5Zw4yM+s4XWlUbaeQtAZwZ/74CdLAodLK0oMj4r1OaZiZmbULB85Wiog3yFNQJJ0MzImI02s9XlKPwrQUMzNrcL5V244kHSZpnKRJks7No3iXlvR2zmQ0BRiclyj7VSHf7QBJt0t6RtLhnf09zMxsEQfOdiKpL2kKynYR0Z/Uuz8o714VuDci+kXEQ3nbzJzv9mHgz6VjyWt3Vqj/iJwycPyiO8NmZtbefKu2/ewGbA2Mz2ltVwBeyPveA64vK39j/ncqsHREzAXmSlooqWeep/ohZw4yM+scDpztR8DFEfHjxTam1VDmxUfnAZUyEy1k8SxFC/HPycysYfhWbfv5B/AlSWtCGn0rydnYzcyWcA6c7SQipgI/A/6RBwHdDqzdHudy5iAzs47jzEFdgDMHmZnVx5mDzMzMOkiXC5yShkqKwgolDUPSyZKOy+9HSpqZ53hOknSMpFGSvlkov42kKZKaXZ/TzMw6RlccrXkwaSmxg4GftrYySUtHxAetblVlx0fEmMK51gYekjSGtL7oOcC3IuL9djq/mZnVqUv1OCX1BHYAvk5ONiBpF0n3SvqbpBmSzs/rYyJpjqQzJU2XdKektfL2sTmzz3jgO5J6Sbor9/7ulLSBpFUlPV+oayVJL0haRtLhOQPQZEnXSlqxlvZHxKvA6aRlxo4EpkREtfVEzcysE3SpwAnsA9waEU8Cb0gamLcPBo4GtgB6A/vl7SsB4yOiD3APi/dQl42IQRFxBnA2cGlE9ANGAWdFxDvAJGDnXH5v4LbcO7wuIrbOmYAeJwXySn5buFX7mbzt/NzO44ETqn3RYuagWbOcOcjMrKN0tcB5MHBVfn9V/gwwLiKezcnUryT1SiElFxid3/+lsJ3CdoBtgSvy+8sL5UYDB+b3BxWO6SvpPklTgeFAnyrtPT4i+ufXVICIWAhcAPw9J5CvKCIuzIF90FprrVWtmJmZtbEu84xT0urArsBnJAXQAwjgb/nfompzcIrb59Zw2huBX+VzDwTuyttHAkMjYrKkEcAuNdRVtDC/zMyswXSlHucw4PKI2DAiekXEJ4GZwI6kFUg+lZ9HHkgaPATp+w/L779c2F7uQRYlaB8O3AeQ88c+AvwBuLmwPNjKwMt5NOzwtvqCZmbW+bpS4DyYjyZOvzZvf4Q0QvVxUjAtlZtLCqrTSL3VU6rUfTTw1ZwB6BDgO4V9o4GvsPit3R8D/wQeAJ5o4fcxM7MG1OUzB0naBTguIvausG9ORPTs+Fa1LWcOMjOrjzMHmZmZdZAlOnBKWpCnckyTdJOk1crLRMTYUm9T0mqSvlXYV7W3KenBNmpjr3wruDSn9OZmyn+YXcjMzBrPEh04Seta9o+IvsCbwFHNlF8N+FYzZQCIiO1a2zgzM+t6lvTAWfQQsF7pg6Tjc/aeKZJ+ljefBvTOvdTfSuqZMwE9KmmqpH0Kx8/J/+6SMwmNkfREziervG+gpHskTZB0m6R1CtsnS5pMlWAuaXVJN+T2PSypX2H3lpIekvSUpMPb9CqZmVmrdInAKakH8DnSvEok7QFsQsoY1B8YKGkn4ETgmdxLPR6YD+wbEQOAIcAZpaBYZivgWFJGn42A7fNUk7OBYRExELgY+GUufwlwdM4cVM3PgIk5G9EPgMsK+/qRRvluC/xE0roVvrMzB5mZdYIlPQHCCpImkXqajwN35O175NfE/LknKZD+q+x4kRIY7ERKOLAeabHpV8rKjYuIFwHy+XoBbwN9gTtyrO1Bmru5GrBaRNybj70c2LNC23cA9geIiLskrSFplbzvrxExD5gn6W7SHwA3FA+OiAuBCyGNqq12gczMrG0t6YFzXkT0z0nUbyPdFj2LFBBPjYgLioUl9So7fjiwFjAwIt6X9BywfIXzvFt4v4B03QRMj4hty87xkQFKLVBrpiMzM+tgXeJWbUT8FzgG+J6kpUlB9Gt5tRQkrSfp48BsUlafklWB13LQHAJsWMdpZwBrSdo2n2MZSX0i4m3gbUmlfLbVMgfdV9qX55q+HhH/yfv2kbS8pDVI6foeqaNdZmbWjpb0HueHImJizuxzcERcLmlz0tqWAHOAr0TEM5IeyNND/g78GrgpJ2MfTx1ZfiLiPUnDgLMkrUq6lr8HpgNfBS7OOXNvLxy2NIt6ryfnMlOA/wKHFcpNAe4G1gR+HhEv1XMtzMys/XT5zEGNRNJ3gPUioupyYS3hzEFmZvVpTeagLtPjbHSS/kwaTPSlzm6LmZm1XMMGTkkLgKmkQTgLgG9HRJPZfDoy96yktYEzgc8CbwHvAb+JiPJE8wBERLXFrM3MbAnSyIODSlmBtgROAk7t7AaV5LmeNwD3RsRGeR7nQcD6Fco27B8nZmZWv0YOnEWrkHp1NJXtp6RamZw39nFJF0maLul2SSvkfRtL+kfO+POopN55e6UMRLsC70XE+aVzRsTzEXF2PmaEpBsl3QXcqeS3Sjl1p0o6MJdbR9K9WpRvd0dJPSSNLJT9bvtdVjMzq1cj94ZKyQ2WB9YhBStYlO3nP5LWBB6WdGMsPsqpYpm8bxPSyNvDJV1NSkLwF2AUcFpEXC9peWCpsgxEAm7MyRL6AI820/4BQL+IeFPS/qQMRluSRso+Iule0uLZt0XEL3P2oxVzufVy/t2q80IlHQEcAbDBBhs00xQzM2srjRw450VEf4A8V/IySX2pLdtPtTIAMyNiUn4/AeglaWVSsLoeICLm5/NWy0C0GEl/JGUCei8its6b74iIN/P7HYArI2IB8Kqke4CtSfMzL87p+26IiEmSngU2knQ28DcWn87yIWcOMjPrHEvErdqIeIjUU1uLxbP99Ade5aPZfpoqUykLUDWlDET982vjiPgzaa7mgEL7jiLlyl2rcOzcGr7XvcBOwL+BkZIOjYi3SD3TscCRwJ+aq8fMzDrOEhE4JW1GygX7BrVl+6krI1BEzAZelDQ0n2+5Qhq/ShmI7gKWl/TNQjUrNnGK+4AD8/PLtUjBcpykDYFXI+IiUoAckG8tLxUR1wI/ohCgzcys8zXyrdrSM05IPb/DImKBpFE0n+2nljLlDgEukHQK8D5wQETcXiUD0Ws5yJ4p6QRgFqmH+f0qdV9PWulkMinv7AkR8Yqkw4DjJb2f6z6UdFv5EkmlP2pOqqHtZmbWQZw5qAtw5iAzs/q0JnPQEnGr1szMrFF0q8Ap6ROSrpL0jKQJkm6R9Ok2qPfY/Ey0uXJjJTX5F46kOa1tj5mZtZ9uEzhztp/rgbER0Ttn+zmJRdNUWpPl51iaHhxkZmZdRLcJnMAQ4P2ybD+TgR6S7ssJEh4DkPQVSeNyRp8LcnICJJ0naXzOOvSzvO0YYF3gbkl35217SHooZyC6pjQqt0jSwTkz0DRJvy7bd2Y+x515FK6ZmTWI7hQ4+5ISHlQyAPhORHw6j6I9ENg+zwFdwKLFqH+YHyb3A3aW1C8izgJeAoZExJA8neRHwG4RMYA0qvf/FU8maV3SWqC7kjIFbV2aCgOsBIyPiD7APcBPKzVY0hE5iI+fNWtW/VfDzMxapJGno3SkcRExM7//HDCQlBYPYAXgtbzvSznV3dKkNIBbkBadLvps3v5APn5Z4KGyMluTbhnPAshTbHYiJY5fCIzO5f4CXFepwc4cZGbWObpT4JwODKuyr5jlR8ClEbHY/ElJnwKOA7aOiLckjeSjGYtKx98REQe3vslAmvdpZmYNojvdqr0LWC73GAGQ1A/YsazcncCwnCEISavnDD+rkALsO0prce5ZOGY2sHJ+/zCwvaSN8/ErVRi5O450q3fN/Pz0YNJtWUg/k1KA/zJwf0u/sJmZtb1uEzjz6in7Arvl6SjTSWt8vlJW7jHSM8rbJU0B7gDWyQOJJpKyEF0BPFA47ELgVkl359uvI4Ar8/EPAZvlcksD70bEy8CJwN2kbEITIuKvucxcYLCkaaRnoKe04WUwM7NWcuagDiJpOeBpoG9EvNOWdTtzkJlZfZw5qMHlpAeTgHPbOmiamVnH6k6DgzpNRIwHNu/sdpiZWeu5x2lmZlYHB04zM7M6OHCamZnVwYHTzMysDg6cZmZmdXDgNDMzq4MTIHQBkmYDMzq7Hc1YE3i9sxvRDLexbbiNbcNtbBvV2rhhRLRo2UbP4+waZrQ0A0ZHkTTebWw9t7FtuI1to7u20bdqzczM6uDAaWZmVgcHzq7hws5uQA3cxrbhNrYNt7FtdMs2enCQmZlZHdzjNDMzq4MDp5mZWR0cOBuMpC9ImiHpaUknVti/nKTRef8/JfUq7Dspb58h6fO11tkgbXxO0lRJkyS1elXulrZR0hqS7pY0R9I5ZccMzG18WtJZktSAbRyb65yUXx/vpDbuLmlCvl4TJO1aOKZRrmNTbWyU6zi40IbJkvattc4GaWOb/l63pp2F/Rvk353jaq3zIyLCrwZ5AT2AZ4CNgGWBycAWZWW+BZyf3x8EjM7vt8jllwM+levpUUudnd3GvO85YM0GuI4rATsARwLnlB0zDvgsIODvwJ4N2MaxwKAGuI5bAevm932BfzfgdWyqjY1yHVcEls7v1wFeI82/b6Tf64ptzJ+fo41+r1vbzsL+McA1wHG11ln+co+zsQwGno6IZyPiPeAqYJ+yMvsAl+b3Y4DP5b/Y9wGuioh3I2Im8HSur5Y6O7uNba3FbYyIuRFxPzC/WFjSOsAqEfFwpN+2y4ChjdTGdtCaNk6MiJfy9unACrkn0EjXsWIbW9GW9mjjfyPig7x9eaA0mrNhfq+baGN7aM3/f5A0FJhJ+nnXU+diHDgby3rAC4XPL+ZtFcvk/1jfAdZo4tha6uzsNkL6Zbs93zI7ohXta20bm6rzxWbq7Ow2llySb439uJW3QduqjfsDj0bEuzTudSy2saQhrqOkbSRNB6YCR+b9jfR7Xa2N0La/161qp6SewPeBn7WgzsU45Z41ih0i4t/5WdIdkp6IiHs7u1FLoOH5Oq4MXAscQurVdQpJfYBfA3t0VhuaU6WNDXMdI+KfQB9JmwOXSvp7Z7SjKZXaGBHzaazf65OBMyNiTisfq7vH2WD+DXyy8Hn9vK1iGUlLA6sCbzRxbC11dnYbiYjSv68B19O6W7itaWNTda7fTJ2d3cbidZwNXEEnXkdJ65N+lodGxDOF8g1zHau0saGuY6FNjwNzyM9ja6izs9vY1r/XrW3nNsBvJD0HHAv8QNK3a6xzcW310NavNnnwvTTwLGngTOkhdZ+yMkex+IPvq/P7Piw+8OZZ0kPvZutsgDauBKycy6wEPAh8oTPaWNg/guYHB+3VSG3Mda6Z3y9Der5zZCf9rFfL5ferUG9DXMdqbWyw6/gpFg202RB4ibTaRyP9XldrY5v+XrfV703efjKLBgfVfS1b/AX8ap8XsBfwJGmU1w/ztlOAL+b3y5NGhD2d/we0UeHYH+bjZlAYqVipzkZqI2k02+T8mt4AbXwOeJP0l/OL5BF2wCBgWq7zHHLmrUZpY/6f0wRgSr6OfyCPWu7oNgI/AuYCkwqvjzfSdazWxga7jofkNkwCHgWGNtrvdbU20g6/1639vSnUcTI5cLbkWjrlnpmZWR38jNPMzKwODpxmZmZ1cOA0MzOrgwOnmZlZHRw4zczM6uDAadbNSJrTwefrJenLHXlOs/bkwGlm7SZnbukFOHBal+HAadZNSdpF0j2S/irpWUmnSRouaVxeQ7F3LjdS0vmSxkt6UtLeefvyki7JZSdKGpK3j5B0o6S7gDuB04Adc8L07+Ye6H2SHs2v7QrtGStpjKQnJI0qrGqxtaQHldZ7HCdpZUk9JP1W0iOSpkj6RqdcSOt2nOTdrHvbEticlIXoWeBPETFY0neAo0k5PSH1GgcDvYG7JW1MSm0WEfEZSZuRVsH4dC4/AOgXEW9K2oWUpaUUcFcEdo+I+ZI2Aa4kZROCtEZmH1LatgeA7SWNA0YDB0bEI5JWAeYBXwfeiYit83JgD0i6PdKSdWbtxoHTrHt7JCJeBpD0DHB73j4VGFIod3VELASekvQssBlpMe2zASLiCUnPA6XAeUdEvFnlnMsA50jqDywoHAMwLiJezO2ZRArY7wAvR8Qj+Vz/yfv3APpJGpaPXRXYhLTeolm7ceA0696K608uLHxeyOL/fyjPzdlcrs65Tez7LvAqqbe7FIsvyF1szwKa/n+UgKMj4rZm2mLWpvyM08xqcYCkpfJzz41ISfrvA4YD5Fu0G+Tt5WYDKxc+r0rqQS4kJQjv0cy5ZwDrSNo6n2vlPOjoNuCbkpYptUHSSi39gma1co/TzGrxL9JKE6uQltiaL+lc4DxJU4EPgBER8W6FRYKnAAskTQZGAucC10o6FLiVpnunRMR7kg4Ezpa0Aun55m7An0i3ch/Ng4hmAUPb4suaNcWro5hZkySNBG6OiDGd3RazRuBbtWZmZnVwj9PMzKwO7nGamZnVwYHTzMysDg6cZmZmdXDgNDMzq4MDp5mZWR3+PywNY/IUZgfXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#feature importance bar chart\n",
    "fis = []\n",
    "for f in range(data.shape[1]):\n",
    "    fis.append(vfi(rf, data, f, 4))\n",
    "fis = np.array(fis)\n",
    "\n",
    "indices = np.argsort(fis)\n",
    "features = X_train.columns\n",
    "plt.title('Variable Feature Importances')\n",
    "plt.barh(range(len(indices)), fis[indices], color='b', align='center')\n",
    "plt.yticks(range(len(indices)), [features[i] for i in indices])\n",
    "plt.xlabel('Importance')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0020059298965339023"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fint(rf, data, 1, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#feature interactions heatmap\n",
    "import itertools\n",
    "\n",
    "data = np.array(X_train.values)\n",
    "fints = []\n",
    "for x in itertools.product(range(data.shape[1]), range(data.shape[1])):\n",
    "    if x[0]!=x[1]:\n",
    "        fints.append(fint(rf, data, x[0], x[1]))\n",
    "    else:\n",
    "        fints.append(0)\n",
    "fints = np.array(fints)\n",
    "\n",
    "rfints = fints.resphape(20,20)\n",
    "plt.imshow(rfints, cmap='hot', interpolation='nearest')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAD4CAYAAABBh0sxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAcfUlEQVR4nO3de5hcVZ3u8e/LJZAQCJJEBJw2OCIqGYikwaAiRA06PjpyiSKiEphjDvooogPKI4rxAiMqgjCCxCOJQGC4BEdRjxCQIEIy0ECuxxtDMIyok2ScSLgFwu/8sXfDTqUrvVbSnV5N3s/z9JOqXW+v2tXV+fWqXXvVTxGBmVlJthnoHTAza+XCZGbFcWEys+K4MJlZcVyYzKw42w30DpRq1BDFmKGJ4R0zBn40Izus78eNtelDKvXxAw89kZ4dMzI9y5gx6dmFDyVH//pM+rBPp0cZ+fKM8PLEXM7vwe4Z2f/OyK5Oj0bGz/Y+WBkRo1u3uzC1MWYodE1IDL8qY+DbM7Kv7ftx1y5LH3LIK9OzJy1Mz172rvQsM6alZ0dPSY7OWZk+7CPpUU74Wkb45MTcgRljnpaRvSoj++P06NqMn+0O8PuetvulnJkVp+jCJOlMSUslLZK0QNLrJJ0qqdfJbWrOzMpTbGGSdAjwTuDAiNgfeCvwMHAqaa+6U3NmVphiCxOwB9WBsacAImIlMBnYE7hN0m0Aki6R1FXPrL5Ybzulh9wRkuZJuk/SdZKGD8SDMrPelVyYbgb+RtJvJV0s6bCIuJDqWOTEiJhY586MiE5gf+AwSfu35iSNAj4HvDUiDgS6gE+13qGkqXWR61qR8e6VmfWtYt+Vi4g1ksYDhwITgWskndFD9L2SplI9lj2A1wCLWjIT6u13SgIYAszr4T6nA9MBOkfIq5vNBkixhQkgItYBc4G5khYDJzRvl7Q31RukB0XEXyTNpOezigTMiYjj+nePzawvFPtSTtK+kvZpbBpHdc7Do8DO9bZdgMeA1ZJ2B/6+kW/m5gNvkPSKeuydJGWcpWNmW1LJM6bhwEWSdgWeAR4ApgLHAT+T9Eh9/Oh+4NdU79jd2fj+6S25KcDVknaob/8c8Nst9FjMLIP8QXE9e7UUMxKzE1KD8PwcLkXOGb+n98OYOfv6QEb2Wx3p2RNT120A70+Pzj8iPTvhgPQsX8/Ipro0I5vxnC2fmZ7tmJKxDxm/YzqFe+s3r9ZT7Es5M9t6uTCZWXFcmMysOC5MZlYcFyYzK44Lk5kVx4XJzIrjwmRmxXFhMrPiuDCZWXFKXis3oHY6ACbckhh+dcbAv8rIjkr/AM45sx9Pyk3K+bD6Cen3v1xp9w/QcWn6MpP5TyVHmZCxFGJCxrirdug9021kRrOJk85Oy12Wcf/MTY927NN75jmfTV9GtFQZy4ja8IzJzIrjwmRmxXFhMrPiuDCZWXFcmMysOC5MZlYcFyYzK86gL0ySjqrbhze/npV0vKTrB3r/zCzfoD/BMiJ+APyg+3rdY+544OqImDVgO2Zmm2zQz5ia6pZMZwEfBDokLam3T5H0Q0lzJf1O0hcGdEfNbKMG/Yypm6TtgauAf4qI5ZLGtEQOBsYCjwP3SPpJRHS1jDGVqkUUHUOAf0i88xszdvSlGdnD0pd5TErtdvM5pd//hD+kZ3lRcnJVzjKTnOUY92Vkd0mPjlyYMe616dHLUg805CwdeUtG9qr06PKMZSZdvUd69UKaMX0ZWBoR17S5fU5ErIqIJ4AbgDe2BiJiekR0RkTn6O37c1fNbGNeEDMmSYcDx7DxjlatUwo31DMr1KCfMUl6ETAD+FBEPLqR6CRJu0kaChzJ+l17zawgL4QZ08nAi4FLpPWOn1zdkrsbmE11lOfK1uNLZlaOQV+YIuKfgX9uc/O5jcv/GRFHboFdMrPNNOhfypnZC8+gnzGliIiZwMwB3g0zS+QZk5kVx4XJzIrjwmRmxdkqjjFtkp2BNydmc1Zu5MhZ6vK2tKUmO92cPuRjGctMOg5IH5fd06M3ZOzv0RldUpjxkfTsOZckR9cmdj4BGDIvMfih9DE3eopxq4wlKR0Zy3JOyFgaNOXEnrd7xmRmxXFhMrPiuDCZWXFcmMysOC5MZlYcFyYzK44Lk5kVx4XJzIrjwmRmxfGZ333hyxnZPdOjqzI+iH9k4hnHj+2cPmbWGe2vyMj+OD169CkZ4378PenZHdPP5p6f0zxhRXqWBxLvP+Os6wnHZNz//RnZkzOy+2Zk2/CMycyK48JkZsUpujBJWle3/F4o6T5Jr9+MseZK6uzL/TOz/lH6MaYnImIcgKS3UX2292EDu0tm1t+KnjG12AX4C4Ck4ZJurWdRiyW9u94+RtKvJH1X0lJJN9ftmp4jaRtJMyV9ZQAeg5klKH3GNFTSAmBHYA+e/4SkJ4GjIuKvkkYB8yX9qL5tH+C4iPiwpGupGmFeWd+2HTALWBIRG7yPtV6L8IwW0mbWt0qfMT0REeMi4lXA24HLVTWPE3COpEXALcBePP/xY8siYkF9+V5gTGO8S2lTlKClRfiwfng0Zpak9ML0nIiYB4wCRgPH1/+Or49B/ZlqVgXQPOtkHevPCu8CJkraETMr1qApTJJeBWwLrAJGAP8VEU9Lmgi8LHGY7wE/Ba6VVPrLWLOtVun/ObuPMUH18u2EiFgnaRZwo6TFQBfw69QBI+KbkkYAV0g6PiKe7fvdNrPNoYgY6H0o0ngpUj8rfsio9HFXrUzPjtw7PZvs3zKyGctnuDojm7EkhZs+nxxdq5y1QekyPls/y4QpicEZGQc8L3o8PXt6ejSnMcbSI9KzY+HeiNjg/MJB81LOzLYeLkxmVhwXJjMrjguTmRXHhcnMiuPCZGbFcWEys+K4MJlZcVyYzKw4LkxmVhwvSWmjs0PRlXjK/tqMTh5DMjqfMCsju09iLqfbxUEZ2Xsysm/uPfKcSzOyV2VkczrATElfFsM5GctiUleP/Dx9SP6Ukf1fGdk/Z2TflB7VZC9JMbNBwoXJzIrjwmRmxXFhMrPiuDCZWXFcmMysOC5MZlacTSpMW6J1t6ST6maWiyQtaTS1nCKp1w99Tc2ZWXk2tRlBv7bulvRS4EzgwIhYLWk4VbsmgCnAEuCRXoZJzZlZYfripVx/tO5+MfAosAYgItZExDJJk4FOYFY9Yxsq6SxJ99Szqumq9JQbL+l2SfdKuknSHn3w2M2sH2zqjKlfW3dL2pbqJPhlkm4FboiIGyPiekkfA06LiC4ASf8SEV+qL18BvLM1J2l74CLg3RGxQtKxwNnASc0HtV6L8D1IngPmdNHY56neM90em5yenZ6YOzB9SI7OyGbJWQ6y4nfp2XNS1+UAn/1IenZOxjKT29Oj3JS2D+/RJclDXjcj/e7nn5ienRAZP6+Z6fvbTl+8lDuEqnX3WJ5v3f0m4FnyWndf2926u+4d93aq1VpvAc6XND4ipvWwLxMlfRoYBuwGLGXDZjP7AmOBOVWHcbYF/tg6UERMp/4/3rmfvIjQbIBsdsPLiJhXz45GA+/g+dbdT0t6iPatu5sv5bpbd58XEU/W4wZwN3C3pDnADGBa877rVt8XA50R8bCkaY37Wy8KLI2IQzbnsZrZlrHZx5j6o3W3pD0lNV91jAN+X19+FNi5vtxdhFbWB8ibL36aud8Ao+vZHZK2l7Rf8oM0sy1qc48xQT+07gbOAL5Rv93/JLCC5z+wYybwHUlPAIcA36V69+1PrP/hG625ycCF9X1sB1xA9bLPzAqzSYUpIrZts30lVRHoydhG7huNy4c3Ln+hke/xU3siYjYwu7Hpc/VXb7kFZH1SjJkNFJ/5bWbFcWEys+K4MJlZcVyYzKw4LkxmVhx3SWmjs3Pb6Orq6VzNHkxObXcBPJAeXbowPZu6GOPH6UNy9N7p2XOXpWc/k7Esh70ysjn+d0b2tRnZjCXjqd11hsQF6YO+7dTk6HtuTh/2uozn7NyMTkBn4C4pZjZIuDCZWXFcmMysOC5MZlYcFyYzK44Lk5kVx4XJzIrjwmRmxXFhMrPiuDCZWXG8JKWNzlcqur6dGM5Y57H2wvTskJiYHh53W1puwXvSx5x5XXr25N4jzzkrI5vRgmbV7N4z3UYek7EPOft7fkb2/sTcLRlj5qw5ylg+M//M9OyEA9KzWuglKWY2SCQVJklHSoq68UBRJE2TdFp9eaakZXWTywWSTpE0S9JHGvnX1W3Htx+4vTazjUn9zO/jgF/W/36hl2yvJG0XEc9s7jhtnB4R1zfua3dgnqTrqTq5/Avw0Yh4up/u38w2U68zprot0huBfwTeV287XNIvJP1E0m8kfUfSNvVtaySdX7cCv1XS6Hr7XEkXSOoCPlG3Df95PXu5VVKHpBGSft8YaydJD9ftlj5ctwJfKGm2pGEpDzAi/gx8A/ga1ZGQRRHxy/wflZltKSkv5d4N/CwifguskjS+3n4w8HHgNcDf8nxH6Z2ArojYj6phcnOGNSQiOiPiPKqW3d+PiP2p2oNfGBGrqbqZdDfnfidwUz27uSEiDoqIA4BfURXKnny98VLu7+pt36n383Tg0wmP2cwGUEphOg741/ryv9bXAe6OiAcjYh1wNdWsCqrW4NfUl69sbKexHao2T1fVl69o5K4Bjq0vv6/xPWMl3VH3rDseaNew8vSIGFd/LQaIiGep2pD/34hY1e6BSpoqqUtS14rV7VJm1t82eoxJ0m5U/d3+TlJQddwN4Cf1v03tzjtobn8sYZ9+BJxT3/d44Of19pnAkRGxUNIU4PCEsZqerb/aiojpwHSoThfIHN/M+khvM6bJwBUR8bKIGBMRfwMsAw4FDpa0d3086Fiqg+PdY3a36n5/Y3uru6iPWVHNgO4AiIg1VB11vwX8uJ6RQdXu+4/1u2nHZzxGMxtkeitMxwE/aNk2u95+D9U7XL+iKlbduceoitYSqtnWl9qM/XHgREmLgA8Cn2jcdg3wAdZ/6fd54N+BO8loPW5mg89GX8pFbHjqcURcWBeT0yLinW2+71M9bDu85frvad8G/HpALdsuAS7pITutcXlKT+PVt82kejloZoVLPY9p6/Mo1XuKKY7rPdLt/IwlKZ+Zn7jMBFib2FFlyCcylpl8Pj3K1zOyOcs2Hk2PjkzsOgLkdT65NiObscxjeeJz1jE94/73TY/OyVhmMmlUevaGjO4+7WxSYYqIucDcNrcN34z9MTPzWjkzK48Lk5kVx4XJzIrjwmRmxXFhMrPiuDCZWXFcmMysOC5MZlYcFyYzK46XpLTT/bmXKV6RPuxncpZNHJ4eHTIvLbf2kIwxd0/PclB6dM6y9OykvTP24U0Z2XdlZP+QkT01PfqS0YnBjCVPOV1Scn62Oc/Z0Q+mZ3l5z5s9YzKz4rgwmVlxXJjMrDguTGZWHBcmMyuOC5OZFceFycyKU+R5TJJGArfWV18CrANW1NcPjoi1A7JjZrZFFFmY6qaU4wAkTQPWRETq6Y5I2rbR9snMBplB91JO0gmS7q5bgF8saRtJ20n6H0kX1B1cDpb0n5LOkbRQ0j2SDpR0s6T/kPThgX4cZtZekTOmdiSNBY4CXh8Rz0iaTtU081pgBPCLiDi1zgIsi4gDJF0EfI+qDflwYCHw3R7GnwpMBeh4MVWD8xSTLkh+DJ84MX3Nwj8lJ6EjsdPekBkZg74/PTp/h/TspIxOMUszlvDsPLn3TLeOKenZLHulR4ekLg/6Xsb9D8vIviQ9OunB9N/x5cpYl9PGoCpMwFupVmV11YVnKPBwfdtaNmzO+aP638XAdhHxGPCYpGclDa+7/j7HLcLNyjDYCpOAyyJivY5nkrYDnoiI1mLyVP3vs43L3dcH22M322oMtmNMtwDvlTQKqnfvJHUM8D6ZWR8bVIUpIhYDXwRuqQ9y3wzkfDiHmQ0Cxb+ciYhpLdevAq7qIbprS+6ljcv/p91tZlaeQTVjMrOtgwuTmRXHhcnMiuPCZGbFcWEys+Jow3MSDaCzc9vo6toxLXzi4+kDz/hIcnS+LknOPpqYm3RE8pAsvTk9u19ON5NPZmTvy8i+Mz26PGf5SqQ/Z/D99Gjq782l6UPekLE0KOPHlWXIwvSsDuDeiOhs3e4Zk5kVx4XJzIrjwmRmxXFhMrPiuDCZWXFcmMysOC5MZlYcFyYzK44Lk5kVx2d+t9HZuUt0dW1wQmobb0wfeNGX07M7p0d5dVps1VO9Z7qNfDDj/o/KyL42I3tYevTcE9Ozn8l4bOe+PGPcjOYJJDaQ4LiMMXMc2D/DLj0gPTsWn/ltZoOEC5OZFWdACpOkkHRe4/ppdcfdjX3PNEl/qBtddn/turHvMbPBaaBmTE8BR3d3O8lwfkSMa3z9T3/snJkNrIEqTM9QNZbc4AMwJI2R9HNJiyTd2lt7JklTJP2bpDmSHpL0MUmfknS/pPmSdqtzcyV9q55pLZF0cP88NDPbXAN5jOnbwPGSRrRsvwj4fkTsD8wCmg2lP9l4GXdbY/tY4GiqLr1nA49HxGuBecCHGrlhETEO+ChwWesOSZoqqUtS14oVazf38ZnZJhqwwhQRfwUuB1rfYD2E59szXcH678U3X8pNbGy/LSIejYgVwGrgxnr7YmBMI3d1fd+/AHZpPUYVEdMjojMiOkePHrIZj87MNsdAvyt3AfCPwE6bOU5r++9ma/Bm77zWk7Z8EpdZgQa0MEXEfwPXUhWnbncB76svHw/c0Yd3eSyApDcCqyNidR+ObWZ9pIROvOcBH2tc/zgwQ9LpwAqgeT7vJyV9oHH9yMz7elLS/cD2wEmbsrNm1v+2miUpkuYCp0VEV0q+c5ii65VpY78j48PXf5rxAf85H0K/dnZabsj1Gfc/PSOb2g0B4Dvp0fkZyxsmTMnYh99kZL+Zkf1URjaxG8DFZ6YP+dHcE3ASrVqZnh2Z8ZxpoZekmNkgUcJLuS0iIg4f6H0wszSeMZlZcVyYzKw4LkxmVhwXJjMrjguTmRXHhcnMiuPCZGbFcWEys+JsNSdYZnsa+ENa9Os54yaOCXBD4jITgKPPTgyenj4m78/IJi6vAODx9OiEjOU+yzOWQnQckp49KSN72d7pWX6XFvtozjKiRzKy30uPjtzgIx37SJvn1zMmMyuOC5OZFceFycyK48JkZsVxYTKz4rgwmVlxXJjMrDh9Vpg2pe13wpiHS1rd0hb8rRvJ7ymp7VkfknaV9NHN2Scz6399OWPa1LbfvbmjpS34Le2CEfFIREzeyFi7UjW7NLOC9WVh2ljb79GSZku6p/56Q719cT2LkaRVkj5Ub79c0qR2dyTpoLqF+I6SdpK0VNLYur34kjqzn6S761nWIkn7AF8F/rbelnXCtpltOX29JOXbwCJJX2vZ/i2qLrq/lNQB3AS8GrgTeAPwe+BB4FCq7ryHAB+havl9qKQFjbGOiYh7JP0I+AowFLgyIpZIGtPInQx8KyJmSRoCbAucAYyt24RvQNJUYCpAxwh6KLE92+++tBzAe07sPdPtuuhID89enpY7K31IDsvI5ix1OTAj++f0aMeK9OwNo9Ozl13Ye+Y5GfvLPYm5nA40Oc/DzhnZzw5Ljp6rjDVHbfRpYYqIv0rqbvv9ROOmtwKvkdR9fRdJw6maWb6JqjBdAkyVtBfwl4h4rM7fERE9rcT6EtVT+yQbthkHmAecKemlwA0R8bvG/bfb/+nUTYs699LW0dfKrED98a5cT22/twEmNI4T7RURa4BfUM2SDgXmUjW4nExa992RwHCqur9j640RcRXwD1QF8qeS3rzJj8jMtqg+L0xt2n7fTNVhFwBJ4+rsw8AoYJ+IeBD4JXAaVcHqzaXA54FZwLmtN0p6OfBgRFwI/BDYn2pSnDOBNbMB0F/nMZ1HVXC6nQJ01geh/x/V8Z9u/w78tr58B7AXVYHqdmjL6QKT64PkT9ezoq8CB/UwI3ovsKQ+PjUWuDwiVgF3Slrig99m5eqzY0wRMbxx+c/AsMb1lcCxbb7vg43Ld9EolhExFxjR5i4vrzPrgNc1to+tt3+Vqmi13l/OpwyZ2QDwmd9mVhwXJjMrjguTmRXHhcnMiuPCZGbFUYRPcO5J5whF14S07Kqb08cdeUzGTtyeHl26Mi33koy7H5nRHYR906Pfn5mePWFexj7kyHlvNqfzSM6Sn9SVGzlLUl6bkX1VenRtxu/CkCnpWc3k3ojobN3uGZOZFceFycyK48JkZsVxYTKz4rgwmVlxXJjMrDguTGZWHBcmMyuOC5OZFceFycyK4yUpbUhaQdUkodUoIHEByKDixzX4vBAe28siYoOeNS5MmSR19bS2Z7Dz4xp8XsiPzS/lzKw4LkxmVhwXpnzTB3oH+okf1+Dzgn1sPsZkZsXxjMnMiuPCZGbF2eoLk6QzJS2tuwQvkPQ6SadKGpbwvUm5Ukg6qqWr8QJJz0o6XtL1A71/qSStq/d9oaT7JL1+M8aaK6nP33LfEvso6SRJi+vf3SWS3l1vnyJpz4Rxk3IDIiK22i/gEGAesEN9fRSwJ/AQMCrh+5NypX4BU6k+WXybgd6XzP1e07j8NuD2zRhrLtA52PYReCnwH8CI+vpwYO+cx9Rfj70vvrb2GdMewMqIeAqea2U+mao43SbpNgBJl0jqqmdWX6y3ndJD7ghJ8+q/kNdJGt7TnZZA0iupPjr/g0CHpCX19imSflj/lf6dpC8M6I72bhfgLwCShku6tf75L27MIMZI+pWk79bP4c2ShjYHkbSNpJmSvjJI9vHFVG0K1gBExJqIWCZpMtAJzKpnbEMlnSXpnnpWNV2VnnLjJd0u6V5JN0naox9+FmkGujIO5BfVX5kFwG+Bi4HD6u0P0ZgJAbvV/25L9Vdm/9Yc1WzrF8BO9fXPAGcN9GNs87i3B7qAY+vrY4Al9eUpwB+BkcBQYAmF/VUF1tXP26+B1cD4evt2wC6N5+MBQPXjewYYV992LfCB+vJcYAJwNXDmYNnH+nfxJmA5MAN4V+O+5zafs+7f3/ryFd3ZZq7+nbgLGF1fPxa4bKCe4+3YikXEGknjgUOBicA1ks7oIfpeSVOpfqn2AF4DLGrJTKi33ykJYAjVy8QSfRlYGhHXtLl9TkSsApB0A/BGqkJWiiciYhyApEOAyyWNpfoPfo6kNwHPAnsBu9ffsywiFtSX76UqBN0uBa6NiLMHyz5GxDpJbwcOAt4CnC9pfERM62FfJkr6NDAM2A1YCtzYktkXGAvMqX9/t6X6AzUgturCBNUTTPWXY66kxcAJzdsl7Q2cBhwUEX+RNBPYsYehRPUf+rj+3ePNI+lw4BjgwI3EWk9uK/Zkt4iYJ2kUMBp4R/3v+Ih4WtJDPP9cPdX4tnVUs8Fud1H95z0vIp4cLPsY1dTmbuBuSXOoZk7TmvctaUeqVwOdEfGwpGm0//1dGhE53QT7zVZ9jEnSvpL2aWwaR/WJAo8CO9fbdgEeA1ZL2h34+0a+mZsPvEHSK+qxd6qP4xRD0ouofnk/FBEba6M4SdJu9TGOI4E7t8gObgJJr6L6674KGAH8V/0ffiLwssRhvgf8FLhWUp//se6PfZS0p6TmH5fu311Y//eyuwitrI95Tm58TzP3G2B0PbtD0vaS9kt+kH1sa58xDQcukrQr1ev7B6jeqToO+JmkRyJioqT7qY4VPMz6/0mnt+SmAFdL2qG+/XNUx69KcTLVQdNL6ul6t6tbcncDs6ne+bkyIkp6GQcwVFL3Sx4BJ9QvbWYBN9Yz3y6q5yxJRHxT0gjgCknHR8SzJe8jcAbwDVVv9z8JrKB6fgFmAt+R9ATVO8/fpTpW+CfgnsaQrbnJwIX1fWwHXED1sm+L85IUW09dXDsj4mMDvS+29dqqX8qZWZk8YzKz4njGZGbFcWEys+K4MJlZcVyYzKw4LkxmVpz/D0Ot1dfMb3naAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.base import clone\n",
    "\n",
    "# Sejong Oh: \"Feature interaction in terms of prediction performance\"\n",
    "def fi_pp(model, X, Y, a, b, mode='classification'):\n",
    "    X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2)\n",
    "    \n",
    "    full_model = clone(model)\n",
    "    full_model.fit(X_train, Y_train)\n",
    "    score = sklearn.metrics.accuracy_score(Y_test, full_model.predict(X_test))\n",
    "    \n",
    "    a_model = clone(model)\n",
    "    a_model.fit(np.delete(X_train, [a], 1), Y_train)\n",
    "    a_score = sklearn.metrics.accuracy_score(Y_test, a_model.predict(np.delete(X_test, [a], 1)))\n",
    "    \n",
    "    b_model = clone(model)\n",
    "    b_model.fit(np.delete(X_train, [b], 1), Y_train)\n",
    "    b_score = sklearn.metrics.accuracy_score(Y_test, b_model.predict(np.delete(X_test, [b], 1)))\n",
    "    \n",
    "    ab_model = clone(model)\n",
    "    ab_model.fit(np.delete(X_train, [a, b], 1), Y_train)\n",
    "    ab_score = sklearn.metrics.accuracy_score(Y_test, ab_model.predict(np.delete(X_test, [a, b], 1)))\n",
    "        \n",
    "    if mode == 'classification':\n",
    "        # err(a)+err(b)-err(a,b)\n",
    "        return (2*score-a_score-b_score)-(score-ab_score)\n",
    "    elif mode == 'regression':\n",
    "        #err(a,b)-(err(a)+err(b))\n",
    "        return (score-ab_score)-(2*score-a_score-b_score)\n",
    "    else:\n",
    "        raise ValueError('unsupported mode '+str(mode))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.00012432012432006534"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inputs = np.array(X)\n",
    "targets = np.array(Y)\n",
    "\n",
    "fi_pp(RandomForestClassifier(max_depth=16, random_state=0, n_estimators=10), inputs, targets, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAD4CAYAAABBh0sxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAcaUlEQVR4nO3de5hcVZnv8e8PghIIASMBSRiISsKMZiCSBkRuiZdRkXNEyYBMUIIc8jgOAnrQQVGMMxp1RkYIo2g4SLhEBgQ8ioKASMvVSToQSDhKYLgMgmK4DBBuA8l7/ti7YafoSq+VdKdWk9/neerpql1vrb12V/Xba+/aa7+KCMzMSrJRpztgZtbKicnMiuPEZGbFcWIys+I4MZlZcYZ1ugOl2nqkYtzotNiV96S3u/GYtetPv7ZNjPt9RpvPZsTu/NqM4OfTQx/PaPaxjNiNM2JXZMSOzYh9buDX/6cn02O3Sfx85/aBt6R/FhYtev6RiHhFT5yY2hg3Gnq+kRb75CHp7Y78u7XrT7+OTYw7MaPN2zJir98xPTaWpcdenNGHCzNit8iIvSEjdnZG7F0Dv/7TrkiPPS7jc5v1O+hJ/yxIy+7va7l35cysOEUnJkknSbpD0u2SFkvaU9LxkjZLeG1SnJmVp9jEJGkv4EBgt4jYBXg38ABwPJCScFLjzKwwxSYmYDvgkYh4HiAiHgGmAWOAayVdCyDpDEk99cjqK/WyY/uI+ytJN0u6RdKPJI3oxEaZWf9KTkxXAX8maZmk70raPyLmAA8BUyNiah13UkR0AbsA+0vapTVO0tbAF4F3R8RuQA/wmdYVSppZJ7me5RnfbpjZwCo2MUXECmAyMBNYDlwoaUYfoYdIugW4FXgr8JY+Yt5eL79R0mLgCOAVXx1ExNyI6IqIrtEjB2Y7zCxf0acLRMRKoBvolrSEKqG8RNIbgROA3SPicUnzgE37aErA1RFx2OD22MwGQrEjJkk7SxrfWDQJuB94ipfPQhkJPA08IWlb4P2N+Gbcb4C9Je1Ut725pAmD2X8zW3slj5hGAKdL2gp4EbibarfuMOAXkh6qjx/dCvyO6hu7Gxuvn9sSNwO4QFLvaalfBDLO9DOz9UW+UFzfut6sSD3zmy9kNPxUeuiyh9NjJ5yVGHh9eps5Z/suujs9dnLq9BlYfQzcj+556bFTDk6PPfOS9NijZ6THsm9i3D4ZbR6fEZsjY8bC5Qemx34AFtVfXq2m2F05M9twOTGZWXGcmMysOE5MZlYcJyYzK44Tk5kVx4nJzIrjxGRmxXFiMrPiODGZWXE8JaWNcVKclBh79K4ZDWdcAP6I1A4A5zyUGJh6AXyA1Ck5UF39KtGTGUUORt6Z0YcJe6THnrYgPfa4nL+RjyVHhs5List5yyZ8LSP4yIzYi9JDuzOmxUz1lBQzGyqcmMysOE5MZlYcJyYzK44Tk5kVx4nJzIrjxGRmxRnyiUnSh+ry4c3bKknTJV3c6f6ZWb6SixEkiYgfAz/ufSxpJjAduCAi5nesY2a21ob8iKmpLsl0MvBRYAdJS+vlMyT9RFK3pLskfbmjHTWzNRryI6ZekjYBfgj874j4T0njWkL2ACYCzwALJf08Inpa2phJVSKKHQRHDx+Ejn4+PfScjIoq7JQYl1HNhDelhz54RXrs2BkZfZhwfkbwsemhOfM8pik9NqOyjFKrn2S0mTM1iO3SS9Bsc3z6G5wzK6adV9OI6R+BOyLiwjbPXx0Rj0bEs8Cl9FEUZ7US4RmfRTMbWK+KEZOkKcDBwG5rCGudienZy2aFGvIjJkmvA84GPhYRa9r5eY+kUZKGAwexetVeMyvIq2HE9AlgG+AMabX9rwta4hYAlwDbA+e3Hl8ys3IM+cQUEV8Hvt7m6W827v8+Ig5aD10ys3U05HflzOzVZ8iPmFJExDxgXoe7YWaJPGIys+I4MZlZcZyYzKw4G8QxprUyAtg7MXbbQerDnIzYxLOyYkx6k/pZeuzYP6bHcmtG7B8OT489KqPdU9NDn9w5PXbkZhl9SJxGNGGLjDYzKtBcrvRpJn/KmWeSMYVmZpsueMRkZsVxYjKz4jgxmVlxnJjMrDhOTGZWHCcmMyuOE5OZFceJycyK48RkZsVRhK8w25edpfh+YuyUjLOIuT4jdnxG7JouKtyUczZ5ztncd+2RHjt+QXps6gX7Ie/3dVF66KKMs6nfkNGFsYln1ncfmN5mzknikzN+B+Sc0X5IeqieYVFEdLUu94jJzIrjxGRmxSk6MUlaWZf8vk3SLZLesQ5tdUt6xZDRzMpT+tUFno2ISQCS3kt1be/9O9slMxtsRY+YWowEHgeQNELSNfUoaomkD9bLx0n6raQzJd0h6aq6XNNLJG0kaZ6kr3ZgG8wsQekjpuGSFgObAtsB76yXPwd8KCKelLQ18BtJP62fGw8cFhFHS7qIqhBmb53pYcB8YGlEvOIKM80S4YN1iSUz61/piam5K7cXcK6kiYCA2ZL2A1YBY3k5l9wbEYvr+4uAcY32vg9c1FdSgqpEODAXqtMFBnhbzCzRkNmVi4ibga2B0cD0+ufkOnE9TDWqAni+8bKVrJ58bwKmStoUMyvWkElMkv4c2Bh4FNgS+FNEvCBpKrBjYjNnAZcDF0kqfbRotsEq/Y+z9xgTVLtvR0TESknzgcskLQF6gN+lNhgR/yJpS+A8SdMjYtXAd9vM1oWnpLTR9TpFz7sSgx/OaHjXjNiz00MvfyYt7oCcqQU58xtypo6kXwMfThqcdpdlXDB/wokZfciZxvP5xLi7MtrMKciQM38mZyrTU+mhOtBTUsxsiHBiMrPiODGZWXGcmMysOE5MZlYcJyYzK44Tk5kVx4nJzIrjxGRmxXFiMrPilD5XrnM2If2U/YsnpLc7e1l67OnpoQd8/PHEyDenN8qe6aHjM+aDjMnoQsa0HG5ND50wIqOqy+yMqi6HpoeyX2JcxueAYzNiM6a6LMu4buyEgzP60IZHTGZWHCcmMyuOE5OZFceJycyK48RkZsVxYjKz4jgxmVlx1ioxrY/S3ZI+XhezvF3S0kZRyxmS+j0TJjXOzMqztidYDmrpbknbU13tebeIeELSCKpyTQAzgKXAQ/00kxpnZoUZiF25wSjdvQ3VJc1XAETEioi4V9I0oAuYX4/Yhks6WdLCelQ1V5W+4iZL+rWkRZKulLTdAGy7mQ2CtaqSImklsIRG6e6IWFTXatusWbqbqmT3jsDdQFdELK5Ld/80Is6X1A2cCBxHXbpb0sZU9d/+ArgGuDQiLqvX3Q2cEBE99eNREfFYff88qkq7lzXjJG0C/Br4YEQsl3Qo8N6I+HjLdr1UInyHEUy+/4jEX8jFGb+81MoYAF/IiE2dipBTxWNaRuwHPpoee9156bH7jUqPnfZYeuyR6aFZ729GhZDk6St/nfE7iPTfwZMZw5KRP0uPzdl30hZ9V0kZiF25AS/dXdeOex+wO/Au4NuSJkfErD76MlXS54DNgFHAHcBlLTE7AxOBqyVBVTjzD60NNUuEd23jEuFmnbLOk3gj4uZ6dDQaOICXS3e/IOk+2pfubu7K9ZbuPiUinqvbDWABsEDS1VTTOWc1112X+v4u1UjsAUmzGutbLRS4IyL2WpdtNbP1Y52PMQ1G6W5JYyTt1nh+EnB/ff8pXi7F2JuEHqkPkDd3PppxdwKj69EdkjaR9NbkjTSz9WptR0yDWrqb6pjTt+qv+58DlgOfqEPnAd+T9CywF3Am1bdvfwQWNppsjZsGzKnXMQw4lWq3z8wKs1aJKSI2brP8Eaok0JeJjbhvNe5Padz/ciP+nW3WcQlwSWPRF+tbf3GLSb8Cjpl1kM/8NrPiODGZWXGcmMysOE5MZlYcJyYzK46rpLRx13J473fSYnNmbhx9fUbwwxmxO2XEpsqZXjE+fZrJmRlTFo4+K32KxZmX9B/zUrupFXAgb/pKTpWS1Onln8+YapNR+WTkrumxHJIROy8jtg2PmMysOE5MZlYcJyYzK44Tk5kVx4nJzIrjxGRmxXFiMrPiODGZWXGcmMysOE5MZlYcT0lp403Ajwaj4dkZsftkxN6aGLfdhPQ2xy9Lj52wsP+Y2tHv3z293YyqLgekh/Jg4nQjgLH/K6Pht2XEPp0YlzHNhCsyYudlxOZUapmdMYWmDY+YzKw4SYlJ0kGSoi48UBRJsySdUN+fJ+neusjlYknHSpov6W8b8XvWZcc36VyvzWxNUnflDgNuqH9+uZ/YfkkaFhEvrms7bXw2Il4qUShpW+BmSRdTVXL5V+CTEfHCIK3fzNZRvyOmuizSPsBRwEfqZVMkXSfp55LulPQ9SRvVz62Q9O26FPg1kkbXy7slnSqpBziuLhv+q3r0co2kHSRtKen+RlubS3qgLrd0dF0K/DZJl0jaLGUDI+Jh4FvAP1FVWrk9Im7I/1WZ2fqSsiv3QeAXEbEMeFTS5Hr5HsCngLcAbwY+XC/fHOiJiLdSleVujrBeExFdEXEKcDpwTkTsAswH5kTEE1TVTHqv2HMgcGU9urk0InaPiF2B31Ilyr78c2NX7i/rZd+r+/lZ4HMJ22xmHZSSmA4D/q2+/2/1Y4AFEXFPRKwELuDl75BWARfW989n9e+WLmzc3wv4YX3/vEbchbxc1f0jjddMlHR9XbNuOtCuYOVnI2JSfVsCEBGrqMqQXxERj7bbUEkzJfVI6mkbZGaDbo3HmCSNoqrv9peSgqribgA/r382tT7ua3nKF6Q/BWbX654M/KpePg84KCJukzQDmJLQVtOq+tZWRMwF5gK8rdpeM+uA/kZM04DzImLHiBgXEX8G3AvsC+wh6Y318aBDqQ6O97bZe7XZv2ksb3UT9TErqhHQ9QARsYKqou5pwM/qERlU5b7/UH+bNj1jG81siOkvMR0G/Lhl2SX18oVU33D9lipZ9cY9TZW0llKNtv6hTdufAo6UdDvwUeC4xnMXAoez+q7fl4B/B24ko/S4mQ09a9yVi4ipfSybUyeTEyLiwDav+0wfy6a0PL6f9mXALwbUsuwM4Iw+Ymc17s/oq736uXkMyGXSzWyweUpKGxtPfi0je3ZMij1C6VM3zsmYMhC3pccqtYrFXel9XZRRpWXykRnTTHLMGZxmx/4sI3jbjNicKSGbJ8btm97kooxKMQ9nVD454P0Z00wuz5hDc9L4PhevVWKKiG6gu81zI9amTTOzXp4rZ2bFcWIys+I4MZlZcZyYzKw4TkxmVhwnJjMrjhOTmRXHicnMiuPEZGbF8ZSUdh5/Hn6UNn3jpJx2t0gP1Z0D3+6yMelNTt4pY/1n75Ee+4MF6bGnpoeO3S89lowqKdySEZtT2Wa3xLiMSjGTD85Y/90ZsX+XEfvzvqeZ5PCIycyK48RkZsVxYjKz4jgxmVlxnJjMrDhOTGZWHCcmMytOkecxSXo9cE398A3ASmB5/XiPiPjvjnTMzNaLIhNTXZRyEoCkWcCKiPhW6uslbdwo+2RmQ8yQ25WTdISkBXUJ8O9K2kjSMEn/JenUuoLLHpJ+L2m2pNskLZS0m6SrJP2HpKM7vR1m1l6RI6Z2JE0EPgS8IyJelDSXqmjmRcCWwHURcXwdC3BvROwq6XTgLKoJAyOA24Az+2h/JjATYDjw4cQqEpfmnK5/cUZsTnWOxHYnZEzx4KiM2GMypplcmR76YMa0ibE523bc+emxRx6eHpszGyN1etLl6U1unlEl5emc6TPtytb25euzMoL7jh1SiQl4N7A70FMnnuHAA/Vz/80ri3P+tP65BBgWEU8DT0taJWlEXfX3Jc0S4a9ziXCzjhlqiUnADyLiS6stlIYBz0ZEazJ5vv65qnG/9/FQ23azDcZQO8b0S+AQSVtD9e2dpB063CczG2BDKjFFxBLgK8Av64PcV5F3JMbMhoDid2ciYlbL4x8CP+wjdKuWuO0b9/9Pu+fMrDxDasRkZhsGJyYzK44Tk5kVx4nJzIrjxGRmxdErz0k0gK7NFD0TEoMzKp+QUXnkmHnpsf+aWFEldk5vUw+lx/JwRmzGFAvOzoidlhH7VEZszjSiIzNiU6d5vD+9ycgo2aOL0mP/PnF6FsA3M07g0cMsioiu1uUeMZlZcZyYzKw4TkxmVhwnJjMrjhOTmRXHicnMiuPEZGbFcWIys+I4MZlZcXzmdxtdmyh6Xp8YnHFR9wdzLhafHsqE1S9b1d7sxDiAXTM6kHOx+t0yYv96VHrsGx5Lj80467l7//TYKavSY/l6YlzGGe3dGWf2T9ksPZaMM8rJmDGg7/jMbzMbIpyYzKw4HUlMkkLSKY3HJ9QVd9f0mlmSHqwLXfbetlrTa8xsaOrUiOl54MO91U4yfDsiJjVu/zUYnTOzzupUYnqRqrDkp1ufkDRO0q8k3S7pmv7KM0maIen/Srpa0n2SjpH0GUm3SvqNpFF1XLek0+qR1lJJewzOppnZuurkMabvANMlbdmy/HTgnIjYBZgPzGk89+nGbty1jeUTgQ9TVen9GvBMRLwNuBn4WCNus4iYBHwS+EFrhyTNlNQjqWd5zrcrZjagOpaYIuJJ4Fzg2Jan9uLl8kznsfqX8c1duamN5ddGxFMRsRx4ArisXr4EGNeIu6Be93XAyNZjVBExNyK6IqJrtL8WMOuYTv/5nQocBWy+ju20lv9ulgZv1s5rPWnLJ3GZFaijiSkiHqM61e2oxuKbgI/U96cD1w/gKg8FkLQP8EREPDGAbZvZACmhEu8pwDGNx58Czpb0WWA5q19F+dOSDm88PihzXc9JuhXYBPj42nTWzAbfBjMlRVI3cEJE9KTEd71Z0fONxMa/ltGRK9JDLx+THntAJFZO+MGy9EY/lR6aVYzglvTQBzOmg4zNKTCQcYH/rGITGdNH4qj+YwDuylj9hFMzgr+QHvrgM+mxf8zoQheekmJmQ0QJu3LrRURM6XQfzCyNR0xmVhwnJjMrjhOTmRXHicnMiuPEZGbFcWIys+I4MZlZcZyYzKw4G8wJltmeI3kuwJm3pTd7dMb0lZwiFoxPnGqSM21jdkZszhSPDGPPygg+MSN224zYnOk2GVPOlfheTJjTf8zarD9nqs3YjKlBW2RMu2rHIyYzK44Tk5kVx4nJzIrjxGRmxXFiMrPiODGZWXGcmMysOAOWmNam7HdCm1MkPdFSFvzda4gfI+niNTy/laRPrkufzGzwDeSIaW3Lfvfn+pay4L9sFxgRD0XEmq66vBVVsUszK9hAJqY1lf0eLekSSQvr29718iX1KEaSHpX0sXr5uZLe025FknavS4hvKmlzSXdImliXF19ax7xV0oJ6lHW7pPHAN4A318v+eQC33cwG0EBPSfkOcLukf2pZfhpVFd0bJO0AXAn8BXAjsDdwP3APsC9Vdd69gL+lKvm9r6TFjbYOjoiFkn4KfBUYDpwfEUsljWvEfQI4LSLmS3oNsDHVpIWJdZnwV5A0E5gJsMPrgTekbfQBaWEALPpOemzGLABILX6SVXIjo0zKxacnh8YN6c3OyIg9J/ZID958QXpsznSbu9ND35E4JeSmnTLWf11G7NkZsRnTTAZgRsrAJqaIeFJSb9nvZxtPvRt4i6TexyMljaCa2bMfVWI6A5gpaSzweEQ8XcdfHxEH9rG6fwAWUs1qay0zDnAzcJKk7YFLI+Kuxvrb9X8u1aiPrnHaMOpamRVoML6V66vs90bA2xvHicZGxAqq/L5vfeumKnA5jbSpiK8HRlBNRdy09cmI+CHwP6kS5OWS3rnWW2Rm69WAJ6Y2Zb+volE+UdKkOvYBYGtgfETcA9wAnEDagPT7wJeA+cA3W5+U9CbgnoiYA/wE2IVqbn1O+UIz64DBOo/pFKqE0+tYoKs+CP3/qI7/9Pp3Xj5Ccj0wlipB9dq35XSBafVB8hfqUdE3gN37GBEdAiytj09NBM6NiEeBGyUt9cFvs3IN2DGmiBjRuP8wjcsJRcQjwKFtXvfRxv2baCTLiOgGtmyzynPrmJXAno3lE+vl36BKWq3r+5t+N8bMOspnfptZcZyYzKw4TkxmVhwnJjMrjhOTmRVHET7BuS9d2yl6ZiQGf/0X6Q1//n3psTnTR54Z+DYXZUyvmPxQemzWFI++zulvY9FR/cf0yik8cs6MjOCzJ6THHpM4jyh9tg98ICP28vS+hlLnPIH2Se+CbmBRRHS1LveIycyK48RkZsVxYjKz4jgxmVlxnJjMrDhOTGZWHCcmMyuOE5OZFceJycyK48RkZsXxlJQ2JC2nKpLQamvgkfXcnfXB2zX0vBq2bceIGN260Ikpk6Sevub2DHXerqHn1bxt3pUzs+I4MZlZcZyY8s3tdAcGibdr6HnVbpuPMZlZcTxiMrPiODGZWXE2+MQk6SRJd9RVghdL2lPS8ZI2S3htUlwpJH2oparxYkmrJE2XdHGn+5dK0sq677dJukXSO9ahrW5JA/6V+/roo6SPS1pSf3aXSvpgvXyGpDEJ7SbFdUREbLA3YC/gZuC19eOtgTHAfcDWCa9Piiv1BswEfg1s1Om+ZPZ7ReP+e4Ffr0Nb3UDXUOsjsD3wH8CW9eMRwBtztmmwtn0gbhv6iGk74JGIeB5eKmU+jSo5XSvpWgBJZ0jqqUdWX6mXHdtH3F9Jurn+D/kjSSP6WmkJJE0ATgY+CuwgaWm9fIakn9T/pe+S9OWOdrR/I4HHASSNkHRN/ftf0hhBjJP0W0ln1u/hVZKGNxuRtJGkeZK+OkT6uA3wFLACICJWRMS9kqYBXcD8esQ2XNLJkhbWo6q5qvQVN1nSryUtknSlpO0G4XeRptOZsZM3qv8yi4FlwHeB/evl99EYCQGj6p8bU/2X2aU1jmq0dR2wef3474GTO72NbbZ7E6AHOLR+PA5YWt+fAfwBeD0wHFhKYf9VgZX1+/Y74Algcr18GDCy8X7cDajevheBSfVzFwGH1/e7gbcDFwAnDZU+1p/FK4H/BM4G/kdj3d3N96z381vfP683thlXfyZuAkbXjw8FftCp93gYG7CIWCFpMrAvMBW4UNKJfYQeImkm1YdqO+AtwO0tMW+vl98oCeA1VLuJJfpH4I6IuLDN81dHxKMAki4F9qFKZKV4NiImAUjaCzhX0kSqP/DZkvYDVgFjgW3r19wbEYvr+4uoEkGv7wMXRcTXhkofI2KlpPcBuwPvAr4taXJEzOqjL1MlfQ7YDBgF3AFc1hKzMzARuLr+/G5M9Q+qIzboxATVG0z1n6Nb0hLgiObzkt4InADsHhGPS5oHbNpHU6L6gz5scHu8biRNAQ4GdltDWOvJbcWe7BYRN0vaGhgNHFD/nBwRL0i6j5ffq+cbL1tJNRrsdRPVH+8pEfHcUOljVEObBcACSVdTjZxmNdctaVOqvYGuiHhA0izaf37viIi91mVbB8oGfYxJ0s6SxjcWTaK6osBTwBb1spHA08ATkrZl9XKNzbjfAHtL2qlue/P6OE4xJL2O6sP7sYh4ag2h75E0qj7GcRBw43rp4FqQ9OdU/90fBbYE/lT/wU8Fdkxs5izgcuAiSQP+z3ow+ihpjKTmP5fezy6s/rnsTUKP1Mc8pzVe04y7Exhdj+6QtImktyZv5ADb0EdMI4DTJW1FtX9/N9U3VYcBv5D0UERMlXQr1bGCB1j9j3RuS9wM4AJJr62f/yLV8atSfILqoOkZ9XC91wUtcQuAS6i++Tk/IkrajQMYLql3l0fAEfWuzXzgsnrk20P1niWJiH+RtCVwnqTpEbGq5D4CJwLfUvV1/3PAcqr3F2Ae8D1Jz1J983wm1bHCPwILG022xk0D5tTrGAacSrXbt955Soqtpk6uXRFxTKf7YhuuDXpXzszK5BGTmRXHIyYzK44Tk5kVx4nJzIrjxGRmxXFiMrPi/H9m5dHJyZEcMQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#feature interactions in terms of prediction performance heatmap\n",
    "import itertools\n",
    "\n",
    "fints = []\n",
    "for x in itertools.product(range(data.shape[1]), range(data.shape[1])):\n",
    "    fints.append(fi_pp(RandomForestClassifier(max_depth=16, random_state=0, n_estimators=10), inputs, targets, x[0], x[1]))\n",
    "fints = np.array(fints)\n",
    "\n",
    "rfints = fints.reshape(20,20)\n",
    "fig, ax = plt.subplots()\n",
    "plt.imshow(rfints, cmap='hot', interpolation='nearest')\n",
    "\n",
    "features = list(X_train)\n",
    "ax.set_xticklabels(features)\n",
    "ax.set_yticklabels(features)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$H^2_{j}=\\sum_{i=1}^n\\left[\\hat{f}(x^{(i)})-PD_j(x_j^{(i)})-PD_{-j}(x_{-j}^{(i)})\\right]^2/\\sum_{i=1}^n\\hat{f}^2(x^{(i)})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#H. Friedman et al.: \"PREDICTIVE LEARNING VIA RULE ENSEMBLES\"\n",
    "from sklearn.inspection import partial_dependence\n",
    "from sklearn.inspection import plot_partial_dependence\n",
    "\n",
    "def h_stat(model, j, data):\n",
    "    jpd, jaxes = partial_dependence(model, data, j)\n",
    "    print(np.array(jaxes).shape)\n",
    "    not_j = []\n",
    "    for i in range(len(data[0])):\n",
    "        if not i == j[0]:\n",
    "            not_j.append(i)   \n",
    "    notj_pd, notj_axes = partial_dependence(model, data, not_j)\n",
    "    # TBC...\n",
    "    \n",
    "h_stat(rf, [0], np.array(X_train.values))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
